
^ ' 
. \ 

L' 



__ - — r^T^^^oTrHrDYiiAMics 

assemblies inc ^^appemdages 

CF OSBITAL A IitH GEOMETRICAL E 22b 

IHIERACTIOMS p HC $12-00 CSC 

l(TR« systems) - — 


^73-20877 

Unclas 

G3/3J 


77?Mf 

SYSTEMS GROUP 


ONE SPACE PARK • REDONDO BEACH, CALIFORNIA 




STUDY OF THE DYNAMICS OF ORBITAL ASSEMBLIES 
INCLUDING INTERACTIONS WITH GEOMETRICAL 
APPENDAGES 


Contract Nos. NAS 8-26131 
NAS 8-28169 


UNIFIED FLEXIBLE SPACECRAFT SIMULATION PROGRAM 

(UFSSP) 


METHODOLOGY REPORT 


December, 1972 

TRW No. 14938-6010-RU-00 


Prepared 



iAAJ 


Ness 

Project Manager 


Approved : 




R. Gluck, Head 
Analytical Dynamics Section 


Approved : 



D. H. Mitchell, Manager 
Dynamics Department 


(Updated to include modifications as of 18 December 1972) 


for 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
Marshall Space Flight Center 
Huntsville, Alabama 35812 


TRW SYSTEMS GROUP 

One Space Park, Redondo Beach, California 90278 



This report was prepared by TRW Systems Group under Contracts 
NAS8-26131 and NAS8-28169 for the George C. Marshall Space Flight Center 
of the National Aeronautics and Space Administration. The work was ad- 
ministered under the technical direction of Mr. Wayne Holland of the 
Astronautics Laboratory and Dr. Gerald Nurre of the Astrionics Laboratory, 
George C. Marshall Space Flight Center. 



Table of Contents 


Page 

I. Introduction 1 

II. Description of the System Model and Notation 2 

2.1 Topological Tree Model Specification 3 

2.2 Definition of Coordinate Systems and Basic Vector 

Quantities . . 4 

III. Synthesis of the Dynamic Equations 11 

3.1 The Sequencing Algorithm 13 

3.2 The Combining Unit 15 

IV. Rigid Combining Algorithm Specification 17 

4.1 The Auxiliary Variables 17 

4.2 The Rigid Combining Algorithm 22 

4.3 Initialization of the Rigid Combining Algorithm 28 

4.4 Specification of the Torque Td^ 29 

V. Flexible Combining Algorithm Specification 31 

5.1 Auxiliary Flexible Quantities 31 

5.2 The Flexible Combining Algorithm 33 

5.3 Initialization of the Flexible Combining Algorithm .... 38 

VI. Disturbance Subroutine 42 

6.1 Prescribed Disturbances 43 

6.2 Rigid-Body Environmental Disturbances 43 

6.3 Flexible-Body Environmental Disturbances 61 

VII. Orbital Subroutine • • 64 

VIII. Control Subroutine 89 

IX. Mass Properties Subroutine 101 

X. Definition of Quantities 105 

Appendix A. Derivation of the Basic Dynamic Equations 115 

A.l. Equations of Motion for a Single Flexible 

Body 115 


ii 



Table of Contents 


Page 

A. 2. Details of the Flexible Comgining Algorithm . . . .131 

A. 3 Details of the Rigid Combining Algorithm 139 

Appendix B. B.l Derivation of Gravity Nominal Forces and 

Torques 151 

B. 2 Derivation of Gravity Generalized Forces 155 

B.3 Derivation of Aerodynamic and Solar Pressure 

Produced Generalized Forces 161 

Appendix C. Modification of UFSSP to Simulate The Dynamics of 

Rapidly Spinning Flexible Systems 166 

Appendix D 184 

References 193 


V 


iii 


I. 


Introduction 


This document presents the complete equations for the Unified Flexible 
Spacecraft Simulation (UFSS) Program developed by TRW Systems for the NASA/MSFC 
under Contract NAS 8-26131. This general purpose simulation program is 
based on an algorithm which utilizes the digital computer to synthesize 
the dynamic and kinematic equations for a topological tree configuration of 
N interconnected bodies (the interconnected system of bodies forms no 
closed loops), the terminal members of which may be flexible. (For illus- 
trative purposes, Figure 2.1 depicts such a spacecraft model (N=15) where 
the possible flexible bodies are shaded.) Necessary input quantities to the 
dynamics subroutine include the mass and inertia properties of each body 
and the flexible characteristics of each terminal member in addition to the 
specification, for each body, of those bodies to which it connects. This 
latter description involves the specification of the number of rotational 
degrees of freedom at each interconnection along with the associated position 
vectors defining these connections relative to the mass centers of the bodies 
involved. These position vectors can be input as time-varying functions 

if desired, thus affording the capability of studying the effects of time- 

✓ 

varying hinge locations. Springs and dampers are assumed to act at each 
interconnection and structural damping in the flexible terminal members 
is included in the form of equivalent viscous damping. 

Figure 1.1 presents the major subroutines of the UFSS program. The 
Dynamics Subroutine is the subject of Sections III, IV and V of this report. 

The Disturbance Subroutine is documented in Section VI; the Orbit Subroutine 
is documented in Section VII; the Control Interface Subroutine is documented 
in Section VIII. 

[Note that although no control laws are implemented at this time, this inter- 
face routine provides for the future addition of any specific or generalized 
control subroutines.] 

Finally, the Modal Subroutine is documeted in Section IX with a list of 
symbols given in Section X. Appendix A details the derivation of the dynamic 
aquations, while .Appendix B contains the derivation of the flexible distur- 
bances . 
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Figure 1.1 Schematic of the Major Subroutines of the UFSSP 
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II. 


Description of the System Model and Notation 

Consider the case of a multibodied flexible spacecraft modeled as a system 
of N bodies interconnected in a topological tree configuration such that 
only the terminal bodies may be flexible. Figure 2.1 below exhibits such a 
configuration where N=15 and the possible flexible bodies are shaded. 



Figure 2.1 Spacecraft Model in a Topological Tree Configuration. 


This figure also exemplifies the method by which the user numerically 
identifies the individual bodies. 


2 



2 . 1 Topological Tree Model Specification 

In particular, whenever a physical system is modeled as in Figure 2.1, 
it is possible to identify a specific body whose attitude relative to some 
external reference coordinate set is of prime interest; this is usually 
that portion of the spacecraft housing the payload and this particular body 
will be denoted as Body 1 of the system. [An example of a possible external 
reference coordinate set for this program is a set located on a user defined 
Kepler orbit such as specified in Section VII.] 

The remaining bodies are labeled in numerical sequence in a fashion 
which denotes the minimum number of interconnections that must be crossed 
in traversinga path from the specific body back to Body 1. The number of 
connecting points crossed is defined to be the level of the body. (Thus, 
there is a unique body having level zero, and this is Body 1.) 

In general, if Body m is connected to Body n and the level of Body m 
is greater than that of Body n, then Body m is defined to be a branch of 
Body n and Body n is defined to be the limb of Body m. A sub-branch of 
Body n is any member of higher level than that of Body n for which the latter 
forms a link in the chain connecting the member to Body 1. 

The rule for numbering the individual bodies in the topological tree 
configuration is simply that: 

min implies that the level of Body m i the level of Body n 

Figure 2.1 clearly illustrates this labeling scheme. There are three 
level one bodies (2, 3, 4), four level two bodies (5, 6, 7, 8), three 
level three bodies (9, 10, 11) and four level four bodies (12, 13, 14, 15). 

From Figure 2.1, the limb-branch relationships given in Table 2.1 
are determined. These relationships, along with the specification of each 
body as flexible or rigid plus the rigid-body and flexible degrees of freedom 
for each body, constitute input to the computer program. 
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It' should be noted in Table 2.1 that Body 13 is entered as a flexible body 
with zero flexible degrees of freedom. This case is allowed to cover the 
option of, anstraining all flexible degrees of freedom for a given flexible 
body for preliminary analyses. 

2.2 

Figure 2.2 depicts a multibodied flexible spacecraft system modeled 
as a topological tree configuration of bodies traveling through space in the 
vicinity of an attracting body. In most cases, motion of the configuration 
is expressed with respect to an orbital reference axis frame although this 
is not mandatory and the reference axis frame may be arbitrarily specified 
if so desired. 
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Several right-handed, orthogonal coordinate frames are used extensively 
in the following development. The pertinent axis sets are described below, 
the inertial reference frame being given with respect to the earth as the 
orbited body. 


x e ; e e (e=l,2,3) — > Coordinates and unit vectors of the inertially 
£ “£ 

fixed coordinate system with origin 0 g at 
the earth's center, e^ and e^ lie in the 
equatorial plane with e ^ normal to this 
plane and pointing northward; e^ is directed 
along the autumnal equinox. 

x r ; e r (a=l,2,3) —+■ Coordinates and unit vectors of the moving 

reference coordinate frame. If the reference 

frame is an orbital one(Section VIII) , then 

its origin 0 lies on the user defined 

Kepler orbit with e^ and e^ lying in the 

r 1 ? 

orbit plane and e^ normal to it. e^ points 
toward the earth's center and e^ forms an 
acute angle with the orbital tangential velocity 
vector. 


R r 


R 


r 


e 


e 

e 

-e 


R 1 



Position vector of the reference coordinate 
center 0 r with respect to the earth's center 

0 . 

e 

Position vector of Body 1 mass center 0^ 

with respect to the reference coordinate center 

0 . 
r 


[Note: repeated subscripts here and throughout this document denote summation 

over the range of the repeated subscript; i.e., 

.1 




= r: 


r 

Si 


+ r; 


-2 


+ r: 


- 3 ^ 
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In what follows, vector quantities will be expressed in tensor (or 
indical notation) form simply by denoting the components. Thus 



In addition, an "outer product" matrix is formed from a given vector 
as follows: 

0 -a4 • a>2 

0 

<j>2 0 

so that, transforming from vector to index notation, 


a0 


g = f 


ag g e • 


Figure 2.3 presents a schematic of a terminal flexible body (Body j) 
and its limb (Body i) . The following coordinate frames and vector quantities 
can now be defined: 

k k 

x ; e (a=l,2,3) — *■ Coordinates and unit vectors of an axis frame 
a -a 

fixed to Body k. If Body k is rigid, then the 
origin of this frame 0^ is located at its 
mass center. If Body k is a terminal flexible 
member, then 0^ is located at the connecting 
point of Body k with its limb. 




(8=l,2,3)->- Coordinates and unit vectors of an axis frame with 
fixed orientation relative to the limb of Body k. 

The origin of this frame is coincident with 0^. 
Normally this frame is used to define some nominal 
orientation of Body k relative to its limb with the 

attitude variables of Body k defining the orientation 

k k n 

of the x relative to the x ° . 
a p 
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Figure 2.3. Two-Body Schematic 
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-r 

U) 


Angular velocity vector of the reference frame. 


-k 

0) 



to 


r 

e 

-a 


k 

to 


a 


k 


e 

-a 




Angular velocity vector of the x^ frame 

Position vector of 0^ with respect to 0^ . 

Position vector from the Body i mass center to 
its connection with Body j . 


— j 1 i 1 

r = r^ e^ »- For Body j a flexible body, r is the position 

vector from the connecting point to an arbitrary 
mass point in Body j when the body is undeformed . 


The various coordinate frames are related through direction cosine 
matrices as follows: 


where 


x j 







a3 3 


j j 0 

= e^ ' 

• .J 

a 3 

-a 

-3 

ji 

= e j ■ 

i 

’ e 

aB 

-a 

-3 

jr 

= e j ' 

r 

a 3 

-a 

-3 


All adjacent bodies are assumed to be connected through a gimbal 
hinge (see Section K 1) with gimbal angular rotations expressed by the 
coordinates 
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ej (01=1,2,3) = 


the gimbal rotations defining the orientation of the 

x^ frame with respect to the x^,° frame, 
a 3 


In addition, once bodies have been combined, the following column 
matrix is used in the sequel: 




for Body j a branch of Body i, this column matrix has a 
rotational component for every flexible and rotational degree 
of freedom of Body i's branches and sub-branches numbered > j 


PJ 


When a terminal body, say Body j, is deformed, the position vector 
= e^ of an arbitrary mass point in the body is given by (see Figure 2.3) 


p-' (r^ , t) = r^ + u j (r^ , t) 

where r* 1 is the position vector of the mass point when Body j is undeformed. 
In this program, u^ is assumed to be representable as a finite sum of 
vector field functions (orthogonal functions) with time-varying coefficients: 

3 


where 


u^ (P, t) = 


V u j 
/ -j a -a 

0=1 


E 

£=1 

hi 


q J t (O (? j ) 


4 (t) X) * 


1 e J 
ia -o 


0=1 


<J>o (r^) = the Jlth orthogonal function describing spatial variation of u^ 

* 


»* 


generalized coordinate describing the time variation of u"* • 


rij = number of terms in the series expansion. 
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III. Synthesis of the Dynamic Equations 

Attention will now be focused upon the detailed description of the 
algorithm utilized by the computer to generate the system dynamic and kine- 
matic equations. The state vector is taken to consist of the scalar elements 






Y 


5 


8 j 

Y 




J 




where 


a, $ = 1,2,3 

j = , • • • ,N 

k = 0,1,..., n^ 

Y = 0,1,. . . »Pj 


with n^ being the number of flexible degrees of freedom of Body j and 
p^ being the number of rigid-body degrees of freedom of Body j with respect 


to its limb (n. 

3 


and Pj may be zero) 


The inductive algorithm used to synthesize the dynamic equations 
for a spacecraft modeled as in Figure 2.1 is based upon an operation which, 
given the dynamic equations for two separate systems, generates the dynamic 
equations when the two systems are coupled together. Specifically, suppose 
the dynamic equations are known for both Systems A and B in Figure 3.1, 
and that System C is a combination of these two through an r (0 < r 5 3) 
degree of rotational freedom interconnection. (Note that only the shaded 
bodies in Systems A, B and C may be flexible.) 
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Figure 3.1. Combining Two Systems to Form a Third 


The algorithm providing the dynamic equations for the combined 
System C is termed the Combining Unit and can be represented schematically 
by a two-input — one-output device as pictured in Figure 3.2. 


System A 

Dynamic Equations 


System B 

Dynamic Equations 


Combining 

Unit 


System C 

Dynamic Equations 


Figure 3.2. Schematic of the Combining Unit 


The actual combining operation is accomplished by eliminating from the 
dynamical equations the forces and torques of constraint between the Systems A 
and B. When applied to an actual spacecraft modeled as in Figure 2.1, 
this Combining Unit is utilized repetitively to synthesize the dynamic 
equations through appropriate interpretations of the Systems A, B and C. 
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Thus, there are two elementary components in the synthesizing 
algorithm: first; the Combining Unit and second, the Sequencing Algorithm 
which specifies the appropriate Systems A, B and C at each application 
of the Combining Unit. 

Let us first consider the latter component. 


3.1 The Sequencing Algorithm 

In order to specify the sequence of combining operations, one need 
merely list for each branch its corresponding limb and the number of 
interconnecting rotational degrees of freedom. Once this information is 
supplied, a computer algorithm determines the sequence of combining operations 
necessary to synthesize the system equations. 

For example, consider the configuration shown in Figure 2.1. From 
the input information supplied in Table 2.1, the Sequencing Algorithm determines 
the step by step procedure outlined in Table 3.1 below, defining the Systems 
A, B and C at every application of the Combining Unit. 

Assuming that the lowest level body in System A has a level greater 
than that of the lowest level body in System B, the development initiates 
at the highest numbered body (15) and considers the interconnection with 
its limb (11) . 


see Section HI. 2 
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Thus, upon the first pass through the Combining Unit, the generated 
equations of motion (C) apply to the two-body system 15-11, where the 
constraint forces and torques between these two bodies have been eliminated 
and the torques along degree of freedom axes are specified by other routines 
of the simulation, such as the control system subroutine. In a similar way, 
the second pass combines body 14 with the system 15—11 to form the three 
body result 15-14-11. This technique is now successively repeated until, 
after the 14th pass, the complete equations of motion (i.e., the second 
derivative of the state vector components) reside in C and are ready for 
the computer’s integration package. This sequence of operations is retained 
by the computer and applied, at each integration step, to form the dynamic 
equations . 


Table 3.1. Synthesis of Dynamic Equations 
for System of Figure 2.1 


Use No. 

System A 

System B 

System C 

Input Parameters 
Required for Body 

1 

15 

11 

15-11 

15,11 

2 

14 

C 

15-14-11 

14 

3 

C 

8 

15-14-11-8 

8 

4 

C 

4 

15-14-11-8-4 

4 

5 

7 

C 

15-14-11-8-7-4 

7 

6 

C 

1 

15-14-11-8-7-4-1 
Store C in Cl 

1 

7 

13 

10 

13-10 

13,10 

8 

12 

c 

13-12-10 

12 

9 

C 

6 

13-12-10-6 

6 

10 

C 

3 

13-12-10-6-3 

3 

11 

C 

Cl 

13- 12-10-6-3-15- 

14- 11-8-7-4-1 
Store C in Cl 


12 

9 

5 

9-5 

9,5 

13 

C 

2 

9-5-2 

2 

14 

C 

Cl 

15 through 1 



C now h< 

juses the sys 


:em dynamic equations 
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Determination of the new translational and rotational positions is 
a considerably simpler task as it is merely necessary to integrate directly 
the rate variables already present in the state vector. This procedure 
applies to all but the body 1 rotations for which the attitude direction 
cosines are desired. These can be obtained by integrating the conventional 
direction cosine equations as will be demonstrated later, [see (4—52) J 

When defining the Combining Unit in the following section, it is 
convenient to identify a system by specifying the lowest leveled body and 
its lowest numbered branch. As shown in Figure 3.3, let Body j of Level 
(N+l) be the lowest leveled body of System A [here, N is arbitrary and 
does not refer to the total number of bodies in the configuration] and let 
Body s of Level (N+2) be the lowest numbered branch of Body j . Let Body i 
of Level N be the lowest leveled body of System B and Body Z of Level (N+l) 
be the lowest numbered branch of Body i. [Note that because the sequencing 
algorithm commences with the highest numbered bodies, it is true that j<&.] 
Following the combining process. Body j becomes the lowest numbered branch 
of Body i and the entire combined system is designated System C. Thus, 


system identification 

is as 

follows : 


System 

A — 


js 

System 

B — 


A 

iZ 

System 

C — 


ij 


3.2 The Combining Unit 

As exemplified in Table 3.1, the Combining Unit is utilized repetitively 
to synthesize the dynamic equations for a combined System C given the dynamic 
equations for its components, Systems A and B. In this program, the dynamic 
equations for the rigid bodies of the system are the standard Euler equations 
while the dynamic equations for the terminal flexible bodies are obtained 
via a Lagrangian approach. Because of dissimilarities in the forms of 
the equations for the flexible and rigid bodies, two distinct combining 


15 





algorithms are utilized in the Combining Unit. Choice of the proper combining 
algorithm depends on whether or not System A is a single flexible body in 
a given combining operation (see Figure 3.4 where the Combining Unit is 
represented by the blocks contained within the dashed lines) . 

The following section contains details of the Rigid Combining Algorithm 
used when System A is not a single flexible body; Section V defines the 
Flexible Combining Algorithm used when System A is a single flexible body. 

Derivations of the governing equations are contained in Appendix C. 

* 

IV. Rigid Combining Algorithm Specification 

Assume that System A in a given combining operation (Use No. in 
Table 3.1) is not a single flexible body. In this case, the Rigid Combining 
Algorithm is utilized to synthesize the equations for System C given those 
for Systems A and B. Prior to the first pass through this algorithm, it 
is necessary to compute certain auxiliary variables that are not elements 
of the dynamic state vector, but that can be algebraically determined 
from this vector. [In addition to the auxiliary variables specified in the 
following section, the gimbal torques and certain quantities associated 
with the flexible bodies should be computed at the same time. These quantities 
are specified in Sections 4.4, 5.1 and 5.3.] 

4.1 The Auxiliary Variables 

Specifically, it is desirable to compute and store for all bodies j 

ir ^ • j j 

the scalar components of A^ , R^ , R^ and . 

Using the notation of Section 2.2 and supposing Body i to be the 
limb of Body j , it follows that 


aB 


Jo 1 


= A JJ0 A Q 

ay y6 


where A^o 1 is the in P ut transformation from the Body i axes to the 

nominal Body j axes. Let us first consider the transformation matrix A^° . 

«jjo can be obtained by three successive Euler rotations 9 (ot = l,2»3)> satis- 

oB “ 

fying the constraints of a gimballed hinge, such that 


this algorithm originally developed in Reference 3. 


(4.1) 
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Figure 3.4. Elements of the Combining Unit 
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( 4 - 2 ) 


A^o = G j3 G ^ 2 
ag ay y6 <$g 


where 





1 0 
0 cos 
0 -sin 


cos 6 


j 


sin 0 


j 


cos 0 


bsin 0 


0 

©^ sin 0^ 
0^ cos 6^ 

0 -sin 0^ 

1 0 

0 cos 0^ 

sin 8^ 0 

cos 0^ 0 

0 1 


( 4 - 3 ) 


( 4 - 4 ) 


( 4 - 5 ) 


Rotational constraints in Equations (4-3) to (4-5) are handled by setting 
the pertinent 3^ identically equal to zero, with the condition that first 
0^ is constrained, than 8^ an( i finally, 0^ if necessary. 

In a similar. fashion, 

Jo 1 = g^° 3 G^° 2 G-j® 1 (4-6) 

ag ay y5 53 
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where the G“|° are identical to the defined by Equations (4-3) 

to (4-5) with 0- 1 replaced by G" 10 (a=l,2,3). The 9^0 are input parameters, 
a . a a 

In most instances the 9^° will be constants; only when generalized relative 

displacement between Body j and Body i is desired will it be necessary to 

input the 0^° as time-dependent quantities and in these instances the 0^° (t) 

♦ . ct ct 

and (t) must be consistently prescribed. 

Angular velocities are determined through the relation 

“i ■ “6° + °«6 5 6 (/ 


where 


(cos 0^ cos 0^) (sin 0^) 0 

(■sin 0^ cos 0^) (cos 0^) 0 


(sin 0^ ) 


j+ j 

G ag = ^aB t^ 036 columns removed which imply no degree of 

freedom; i.e., if 0^ 5 0, then the ath column of 

i "1 

G^ is deleted in forming G^ , 

. o 

• j • 2 

while 0^ is equal to 0^ with those rows deleted which imply no degree 

of freedom. Thus, if 0^ has two degrees of freedom, then m* 1 is given bv 

p a 

Ceos 0^ cos 0^) (sin 0^)1 ( 0^ i 


“a = A ag° W B° + ^ sin e 3 cos ) (coS ®3 ) ®2 ' 


(sin 0^ ) 


5 

(4-7) is not coded. The proper expression is (4-10). 
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In case the 0^° (ot=l,2,3) are constant , Equation (4-7) becomes 


or 1 = A J J 
a ag 


"e + 


ej 

ag 6 


In the general case where the 0^° are time-dependent input functions, 
Equation (4-7) becomes 

J 


= A 


Ji 


a). 


+ 


,j+ u 


A lJo qJo 


ag “0 ag “g ' “ay ~yg 

j 


H° 


where G J ° is defined similarly to G „ with 0^° replacing 0^ . (There 
ag . , ctg a ^ 6 a ' 


is no need to define a G 


jo+ 

ctg 


since it is sufficient to set 


Aj 0 . 


for any component 0^° which is constant.) 


It is evident that A^g can be computed directly from the state 
vector and the input parameters 0^° (a=l,2,3). Adopting a procedure 

whereby the auxiliary variables associated with bodies of lower level 
are computed first, it can be assumed that A^ is already known. A^ 
then follows simply from the relation 


A jr = A jl A ir 
ag ay yg * 

To complete the description of the auxiliary variables. 


y 


A J ! ' A Jr A r * 
ag ay yg 


= R 


R 1 

a 


+ A 


- A' 


> ir (T) 
ag 

Jr(T) 

ag 


. ir (T) 
ag 

*6 

- a^ ( t) 

ag 

B 

.ir (T) 
ag 

£ ij 

B 

- A jr(T) 
ag 

i ji 

g 

~i 

“BY ' 

~r 

U) 0 

ag 

A lr(T) f * ij 

by ( Y 


By 

~r 

b) * 

ag 

A^< T > l^ 1 

By | y 



(4-9) 

(4-10) 


(4-11) 

(4-11.1) 

(4-12) 

(4-13) 
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where the super-dot above a variable denotes time differentiation in the 
designated axis frame. In the above fashion, the auxiliary dynamic variables 
for all bodies can be obtained. It is merely necessary to consider them 
in the proper numerical sequence beginning with Body 1. 

4. 2 The Rigid Combining Algorithm 

Suppose that System A has Body j as its member of lowest level, 

System B has Body i as its member of lowest level (i<j) and that A is 
to be connected to B to yield C as in Figure 3.3. Referring to this figure, 
it is also assumed that Body £ is the lowest numbered branch of Body i in 
System B, s is the lowest numbered branch of Body j, and j<£ . 

The dynamic equations for System A are as follows: 


System A 


All 

’k£ 

A12 

B kg 

_ A13 

\g 

A21 

a£ 

b A22 

ag 

tj A.23 
B ag 

A31 

a£ 

r A32 

B cB 

.A33 

e «b 


Here, 


b A22 b A23 b A32 


ag ag ag 

a and^ g running from 1 to 3; 




"1 

r b 


. n A33 
and B „ 

• o3 


,Al 

- k 

„A2 

"a 

,A3 


(4-14) 


are 3x3 matrices with the indices 


„A2 

■J 

a 


and 


„A3 

"a 


are 3x1 matrices; 


B 


A12 


•j s All 

if O': has M components, then B, is M x M, 

^ a?i A31 A1 

are M x 3, B and B are 3 x M while C. is M x 1 
ai a£ k 


kg 


and B 


A13 

kg 


System B has a similar description and interpretation. 


System B 


b b11 

mn 

b B12 

mg 

B B i 3 

mg 


r~ 

e ,lA 

n 


"c B1 

m 

b b21 

an 

B22 

B ag 

B23 

B ag 


“B 

= 

c B2 

a 

b B3! 

an 

b B32 

B ag 

r B33 

B ag_ 




c B3 

a 


(4-15) 
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Finally, the output of the Rigid Combining Algorithm has a similar appearance. 


where 


System C 


B C11 

pq 

» cl2 

V 

_ C13 

V 


/V 

q 


C C1 

p 

B C21 

otq 

C22 

R C23 


•i 

= 

c C2 

a 

B C31 

aq 

r C32 

R C33 

B cte 

— 


— — 


c c3 

a 



with the column matrix 0^ including only the degree-of-freedom components 
at the interconnection of Bodies i and j . 


To complete the description of the Rigid Combining Algorithm, it 
is now necessary only to explicitly specify the elements of the matrices 

c c 

B and C and to show how this inductive process can be initiated at 
the bodies of highest level by direct use of the computer input data. 
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where 
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where the superscript 0 is again applied as before, meaning that only 

1i 

those components of T along degree-of-freedom axes are pertinent 
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This completes the definition of the Rigid Combining Algorithm. 
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The kinematic equations for Body i (ii2), that is, a determination 
at time t of its attitude rates, is already known since terms such as 0 3 
are already elements of the dynamic state vector. For Body 1, a direction 
cosine approach is desired and the kinematic equations are as follows: 
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•lr 1 .lr , »lr r 

ag ay yg ay yg 

The above relations were coded in the initial UFSS Program. However, 

llT ® 

a much more efficient calculation of A . is realized by utilizing an 

Ctp 

"Euler parameter" technique as follows. 

Utilizing (4-52) , it follows that 

;el .el ~1 
A . = A m . . 

ag ay yg 

©1 

The initial orientation for A . is obtained from 

ag 

A e J: (t ) = A er (t ) A r J(t ) < 
ag v o ' cty o yg v o' 

The matrix A „ can be written in terms of A _(t ), a constant, and a 
ag ag ° 

matrix C _ , a function of time as follows 
ag 


A 6 ^(t) = A e± (t ) C Q (t) , 

ag ay o yg v 


.el 


(4-52-a) 


el 

so that C Q (t)- also satisfies the same differential equation as A ; i.e., 
Yp otp 


C Q = C 0) „ 
ag ay yg 


Instead of solving the above matrix differential equation, C 
can be expressed in terms of the four Euler parameters consisting of a scaler 
X and a vector resulting in four scalar differential equations to be 

integrated. 

To introduce the Euler parameters, note that any orientation of 
a body may be achieved by a counterclockwise rotation about an appropriate 
unit vector through an angle 6 . Accordingly, has the repre- 

sentation 


C „ = <$ „ + sin 9 e n + (1 - cos 0) a e . 

ag ag ag ay yg 
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The Euler parameters are defined in terms of e^ and 0 by 


X = cos (0/2) 


k = sin (0/2) e 
a a 


6 q + 2x^o + 2k k. 
ag A ag ay yg 


The differential equation for C^g leads to corresponding differential 
equations for X and , namely 

; - - 1/2 » im K 

a a 


■ 1/2 x (“i - “It *s) ■ 


2 (T) 

By definition, x + K < is equal to one and indeed this function is 
./ a. a 

an integral of (4-52-c). This fact is used to provide a check on the com- 
putation through calculation of 

4 - |1 - X 2 + * <T> K | . 

1 A o a 1 

As seen from (4-52-a) , the initial value of C^g is the identity matrix. 
Thus, from (4-52-b), initial conditions for (4-52-c) are 


x(t )= 1 ; k (t ) = 0 . 

o a o 

Finally, the equations (4-52-c) and (4-52-e) are solved for x(t) and K (£)- 
with (4-52-b) subsequently used to obtain C a g(-0 • From (4-52-a) 

A^-gCtX = C (T) (t) A X ®(t ) . 

ag ay yS o 


Thus, since A^ = A^ r A r ® , A^ is obtained from 
ag ay yg ag 
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(4-52-f ) 


A^CtX = C (T) Ct) A^(t ) A*®(t ) A^g (T) Ct). 
otg ay ycr o ap o pg 


It now remains to detail the initialization of the inductive process. 

4. 3 Initialization of the Rigid Combining Algorithm 

When one of the inputs to the Rigid Combining Algorithm is a single 

A A 

rigid body (say Body k is the System A) , then the matrices B and C 
are initialized as follows: 

,A22 T k 


B 


ag 


ag 


D A33 k .kr 

B . = m A . 

ctg ctg 


- “ A f - o 

ctg ag 


All the remaining sub-matrices of B are void. In addition, 


c A1 

ot 

is 

void 




c A2 

ot 


~k 

0) 

ay 

I* 

y6 

“5 + 

_ke 

T 

ot 

c A3 

a 

= 

k 

- m 

A kr 

ay 

[ 2 V 

R k 

R 6 


, , T r ~r ~r s _k 

+( “yj + “ye V *5 


+ F 


ke 


If Body k is, instead, a single rigid body treated as a limb (System B) , 

A A 

then the above specified sub-matrices of B and C define the sub- 
B B 

matrices of B and C . In either event, the Rigid Combining Algorithm 

can be readily initialized in its inductive procedure of synthesizing the 

1&& k© 

system dynamic equations. For specification of and , see Section 6. 
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(needed in (4-46)) 


4.4 Specification of the Torque T^ 

The relative rotations 0^ of Body j with respect to its limb, 

Body i, are assumed to take place about a gimballed hinge nominally aligned 
with the axis frame (x^°, x^ 0 , x^ 0 ). The components of the torque T^ 
are functions of the reaction torques at these gimbal axes. 

Specifically, let T^ , T^ , T«^ be the gimbal hinge reaction 

gag 

torques about the x® , x® , x® gimbal axes of Body j . These three 
gimbal torques may be arranged in a 3 x 1 column matrix such that 



where the elements of T^ 1 are not orthogonal components of a resultant 
torque vector. 

It can readily be shown that the torque components of the total 
gimbal reaction torque, transmitted from Body i to Body j, defined in 
Body j coordinates are given by 
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(4-60) 


29 



(It is assumed that the gimbal rotations are such that 0 £ 9^ < u/2 

so that cos 9^ 4 0.) 

The gimbal torques about the gimbal axes of Body j are 


T jh = + T jd + T jm 

a a a a 


where 


pJ s = 


- K jl 0 j 

K a0 8 a 


J2 

e e 


k 33 e 33 
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are the spring restraint torques. 


r-5 d = 


r n -j 

~ C aB e B 


r j2 iJ2 
C aB 9 e 


r j 3 Aj3 

°ag 6 B 


are the damping torques and 


(4-61) 


T^ m is the array containing the motor torques about the 
gimbal axes of Body j in case Body j is a controlled 
body. 

In the above relations. 


K 


Jn 

«B 


kj n 0 0 

0 k^ n 0 

0 0 k^ n 


where k| n , k^ 11 , k^ n are the spring constants associated with the x^ , 
x® » x 2 gimbal axes of Body j in units of FL/ (rad) n ; 
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J n 


0 


Jn 


J n 

'2 

0 


0 

J n l 


• , g 

where c Jn , c 3n , c 3n are the damping constants associated with the x x , 
1 ~ 2 3 . _ , ,,^n ...n.i, 


*2 ’ 


x® gimbal axes of Body j in units of FL/(rad/T) . In addition, 


the following notational convenience is utilized: 


e 3 . 11 


(e 3 ) n 


(e 3 ) n 


(e j 3 ) n 


The body j control motor torques T 3 ™ are computed in the control 
routines . 

V. Flexible Combining Algorithm Specification 

Consider now the case where System A in a given combining operation 
is a single flexible body. In this case, the Flexible Combining Algorithm 
is utilized to synthesize the dyriamic equations for System C given those 
for Systems A and B as shown in Figure 5 . 1 . Prior to the first pass through 
this algorithm, the following quantities are computed. 

5.1 Auxiliary Flexible Quantities 

The following quantities are computed for each flexible body 
at the time the auxiliary variables are calculated: 
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(5-4) 
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In all the above quantities, 

a, 8 = 1.2,3 


k, £ 1,2,.. . ,n . 
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(f 

(f 

a 

0 

(‘ 


The arrays 2“^ , an< ^ ^la8 are °^ ta ^ ne< ^ fro® the Mass Properties 

Subroutine and will be defined in Section IX. 


5. 2 The Flexible Combining Algorithm 


The Flexible Combining Algorithm is used only when System A in a 
given combining operation is a single flexible body, call it Body j. Once 

again, assume that System B has Body i as its member of lowest level (i<j) 

and that System A is to be connected to System B to yield System C as in 
Figure 5.1. Referring to this figure, it is also assumed that Body l is the 

lowest numbered branch of Body i in System B and that j<£. . 


-5) 
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- 8 ) 
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- 10 ) 

i-ll) 
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The dynamic equations for System B will be specified first since 
they are identical to the System B equations for the Rigid Combining Algorithm. 


Once again, 
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matrices; if 
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(5-13) 


Specification of the System A equations will now be given. It should 
be noted that since System A is a single body, the following form of System A 
occurs in the synthesizing algorithm only as a matrix loading operation 
and never as the result of a combining operation. 


Here, 
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Thus , B 
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(Note here that the sub-matrices occurring in the equation for a PP ear 


with a superscript 
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and B 


denoting removal of the Xth row of the sub-matrices 
0 ; hence, 5 and y need not run fully from 1 to 3 as must a 
However, the full sub-matrices must be loaded as they appear 


in the System C equations except for Equations 5.17-5.20) 

Finally, the output of the Flexible Combining Algorithm has a form 
identical to the output of the Rigid Combining Algorithm. 
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"i 

again with the column matrix 0^ including only the degree-of-freedom 
components at the interconnection of Bodies i and j . 

To complete the description of the Flexible Combining Algorithm, 

it is now necessary to explicitly specify the elements of the matrices 

Q Q A A 

B and C and to show how the elements of the matrices B and C 

are initialized. (Note that if System B is a single rigid body, it is 

initialized as shown in Section 4.3.) 
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( 5 - 21 ) 
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This completes the definition of the Flexible Combining Algorithm. 


It now remains only to detail the initialization of Systems A and B for 
this combining algorithm. However, before so doing^ the expanded form of 
System C will be presented for completeness. 
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5.3 Initialization of the Flexible Combining Algorithm 

As stated previously, if System B is a single rigid body, then its 
initialization is identical to that of the Rigid Combining Algorithm presented 
in Section 4.3. Initialization or loading of System A is as follows: 
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where and are respectively the Body j generalized mass and 

stiffness matrices obtained from the Mass Properties Subroutine while 

is the input generalized modal damping matrix. 
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In the above equations, , Qy ^ and & are externa l gen- 

eralized forces to be specified in Section VII. Finally, the through 

S jl ° are defined as follows (these quantities are computed for all flexible 
bodies at the time the auxiliary variables are calculated) . 


i^ 1 

a 

= 

A jr ( 

• i 

2 % + 

•r 

“aB 


+ 

A ji 

aa 

? i 1 . 

ap pB 6 


;3 2 

a 

= 

£ ij 

a 

+ 2 ; a 6 


;^ 3 

a 

= 

A jjo 
a a 

g j ? +/a j1 
aB \ aa 

~i 

to 

ap 


+ 

A jJo 

aa 

ej° = P 

aB 8 

,100 

a 

J4 

’kB 

= 

yj 

kB 

- H jl 
H kB 


. j 5 
* 

a 

= 

- 

aa 

J 1 1 4 

A aS “8 + 

G j+ 

ay 

46 

'a 

= 

d^ + 
a 

3 

a 


; J 7 

'a 

= 

s^ 1 

a 

+ A ji S j2 
aB S 8 
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'al 

= 

yj(T) 

aSL 

4- H^< T) h 

a£ 

a i 

49 

aB 

= 

1^ ^ 
aB 

+ + hJJ° 

aB aB 

410 

k 

= 


+ “&> *4 

4 


3 i ~ r ~ r D i| 

l. + to to . R. I 

8 ap pB By 


- -j 

pe 

[see (4-50)] 


to- 

aa at 


io \ pio I A jo 

t ) e * e J 9 e 


•j° 

0 J 


+ S* 3 HP 17 


(5-49) 


(5-50) 

(5-51) 

(5-52) 

(5-53) 

(5-54) 

(5-55) 

(5-56) 

(5-57) 

(5-58) 

(5-59) 
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In addition, the following quantities will be needed for the dis- 
turbance calculations to follow and should be calculated with the above; 


ka 

_ r(T) je(T) 

a e 

B ka<$ 

(5-60) 

ojl2 

k£B 

= r (T) je(T) 

a S A <Sa 

kilaB 

(5-61) 

oj 1 3 
S £a 

= q j S jl2 
^k k£a 


(5-62) 

ka 

= B j A Je a r 

kag Bp p 


(5-63) 

ag 

= H^ 2 < T > + H *? 

ag aB 


(5-64) 


Once again the arrays and Cjj^ are obtained from the Mass Properties 

Subroutine, while a^ is determined by the Orbit Subroutine. 
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VI. Distubance Subroutine 


The Disturbance and Control Subroutines supply the perturbing forces 
and torques required in Equations (4-56), (4-57), (5-39), (5-44) and (5-49). 


At present, the Control Subroutine is not implemented, although 
aralized Control Interface Routine described in Section VIII is 
part of the program. 

the Gen- 
an integral 


The general disturbance 

equations are as follows: 



a 

_ T jD 

a 

+ 

T JC 

a 


(6-1) 

F- 3 ® 

a 

= 

a 

+ 

F^ C 

a 


(6-2) 

of 

- 

+ 

of 


(6-3) 

a 

- 

01 

+ 

r! c 

a 


(6-4) 

®f 

“®f 

+ 

®f 


(6-5) 


where all the superscript J quantities are presently void control vector 
forces and torques and 


T j° = 

x jA 

+ 

T- 3 S + 

+ 

t J g + 


(6-6) 

a 

a 


a 

a 

a 

a 

F j° = 

f Ja 

+ 

F- 3 S + 

F- 3 S + 

F^ g + 

F^ P 

(6-7) 

a 

a 


a 

a 

a 

a 


of - 

of 

+ 

+ 

CO 

O’ 

of 



(6-8) 

r^ d - 


+ 

R^ S + 

Ri° 



(6-9) 

a 

a 


a 

a 




®f- 

® JA 

+ 

®f 

+ ®f 



(6-10) 
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Definitions of the individual terms in the above relations are amply 
provided by the flow-diagram of Figure 6.1 as well as by the table of 
Section X. 

The following sub-sections present the equations necessary for program- 
ming the disturbances listed above (other than control disturbances). 

Section 6.1 details initial expressions for the prescribed disturbances. 
Section 6.2 (pages 45 through 60) details the environmental disturbance 
equations for rigid bodies, while Section 6.3 presents the environmental 
disturbance equations for flexible bodies. Since the rigid-body environmental 
disturbances are identical to those found in the TRW GSS program, their 
derivation can be found in Reference 1 and is not included in this document. 
However, the derivation of the flexible-body environmental disturbances is 
presented in Appendix B. 

6.1 Prescribed Disturbances 


Prescribed forces and torques are allowed to act about the Body j mass 
center. Specifically, the following analytical representation is utilized 
initially: (t^_ is a reference time, normally equal to zero) 


F jP = C jfl 
a a 

+ c jf2 

a 

(t - t ) + 

r 

c j f3 

a 

sin 


(t - t ) 
r 

(6-11) 

T i p = c j t;L 

a a 

+ c ^ 2 

a 

(t - t ) + 

r 

c jt3 

a 

sin 

u> jt 

(t - t ) . 
r 

(6-12) 

In addition 

to the above representation, 

F^ P 

a 

and 

i P 

T J can be input in 

a 



tabular form. 

6. 2 Rigid-Body Environmental Disturbances 


This sub-section details the equations for the rigid-body environmental 
disturbances. The format herein differs from that of the body of the report 
in that an original document is included here in its entirety. The reason 
for this inclusion is simply that the coding of the program proceeded from the 
equations contained herein and therefore it is felt that its inclusion would be 
preferable to a complete reordering of the equations and updating of the notation. 
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1. INTRODUCTION 


In what follows, forces and torques ^ \ due to environmental disturbances 
resulting from aerodynamic ^ ^ and solar radiation pressure, gravitational 

attraction, and magnetic interaction, acting upon the rigid bodies of the 
spacecraft model are presented. 

To this end and in view of the fact that the magnitudes of the forces 
and torques resulting from aerodynamic and solar radiation pressure are 
dependent on the shape of the surfaces over which they act, four basic 
shapes are provided for use in their computation. These consist of: a 

sphere, flat plate, circular cylinder, and a rectangular parallelepiped. 

In contrast, the gravitational interactions are dependent on the mass 
distribution of the bodies and their relative positions, whereas the mag- 
netically induced forces and torques are only dependent on the positions and 
orientations of the bodies. Therefore, the latter two disturbances can be 
applied to an arbitrarily shaped body. 

In the next section the force and torque (moment) equations for the 
four basic allowable shapes are detailed. The equations are expressed in 
terms of a unit flow vector 6^ , a force P and appropriate constants 
h*^ . For solar disturbances, <5^ becomes a unit vector 6 

" f 

directed along a line from the sun toward the earth, P becomes the solar 
radiation pressure P , and G^ , become functions of the solar 

reflectivity coefficient • For aerodynamic disturbances, 6 becomes 

a unit vector 6 a directed along the negative of the spacecraft's velocity 
vector, P^ becomes twice the dynamic pressure P 3 , and , H* 1 

become functions of the so-called aerodynamic reflection coefficient o . 
More specific values of P , 6 » G^ and are given in Section 3. 

These forces and torques will contribute to the vectors F^ and T^ , 
defining the externally applied forces and torques, of Section IV. 

^ ^ The aerodynamic interaction implied here is that corresponding to free 

molecular flow. 
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Also included in this section are the gravitational and magnetically induced 
forces and torques. Section 4 contains a detailed input/output specification. 

In this appendix, matrix notation is used instead of the index notation 
employed throughout the body of the report. Given the following notational 
clarifications, this use of matrix form should not present difficulties in 
interpretation: 



denotes the coordinate system in which the components of f are 
given: i. e . , 



j equal to r implies components in the orbital reference coordinate 
r it r 

system (x^ , x^ , x^) , while j = e implies components in the inertial 
reference system (x^ , > x ^) • 

(ii) A ~ above a symbol denotes an outer product matrix: i.e.. 



(iii)A superscript T on a matrix denotes its transpose. 


46 



2. GENERAL FORCE AND MOMENT EQUATIONS 


2.1 Forces and Moments on a Sphere 

If Body j is modeled as a sphere, then its surface is defined by 

1 g 

the- scalar radius R J . 

The total force on the sphere is given by 

{f^ - * P f [f + cJ £ l {<5^ (1) 

and the total moment about the sphere's mass center is 

n-nn 

Where L^ P is the vector position of the geometric center of the sphere 
as measured from its mass center 

2.2 Forces and Moments on a Flat Plate 

If Body j is modeled as a flat plate, then its surface is defined 
by its scalar area A^ and by a unit normal vector as shown 
in Figure 1. The plate need not be rectangular. 



Figure 1. Surface Specification for a Flat Plate 
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(3) 


The total force on the plate is given by 

- P f A j | cos n j | [h jf cos v? {e 3 ^ + G jf {^jj] 

where 

“* " j ■ { e3f }j {‘ f }j • (4> 

The total moment about the plate's mass center is 

<5) 

where L 3 ^* is the vector position of the geometric center of the plate as 

measured from its mass center. 

2.3 Forces and Moments on a Circular Cylinder 

If Body j is modeled as a right circular cylinder, then its surface 

i c 1c 

is defined by its scalar radius and length, R J and h J respectively 
and by a unit vector directed along its axis, f as shown in 

Figure 2. 



Figure 2. Surface Specification for a Cylinder 
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The total force on the cylinder is given by 

i F J f l = p f r^ c [*h^ c sin (4/3 H js + 2G js ) + * R Jc (cos | 

1 L (6) 

+ P f R^ C cos 4/3 H^ S h^ c sin " R^ C ' |cos | H^ e J |e^ ^ 

and the total moment about the cylinder's mass center is given by 
{ xJf }j “ I hjC ?f cos ^ (GjS " gJ6) }j 

(7) 

where is the vector position of the geometric center of the 
cylinder as measured from its mass center. 

In the above expressions, G^ S and H^ S pertain to the sides of the 
cylinder while G J and H J pertain to the ends of the cylinder and 

cos n J = { e ^ f }j T 

( 8 ) 

sin = + ^1 1 - cos 2 n^ 

2.4 Forces and Moments on a Rectangular Parallelepiped 

If Body j is modeled as a rectangular parallelepiped, then its 

112 113 123 

surface is defined by three unit vectors (e J , e~ , er ) normal to 
any three nonparallel sides, and the lengths of these three sides 
L J , L J . Orientation of the unit vectors is shown in Figure 3 y 
where in general is normal to those sides with area L^ a L^. 
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Figure 3. Surface Specification for a Parallelepiped 


The total force on the parallelepiped is given by 


= P f I/ 51 L J2 

|cos n 312 | 

j^ 12 cos 

r, 312 {e 312 }. 

,- 3i2 ni 

+ P f L^ 1 l/> 3 

| cos n 313 | 

£ H^ 13 cos 

n 313 {e 313 }. 

I + gJ13 H] 

+ P f L 32 L^ 3 

| cos n^ 23 | 

[V 23 cos 

3 23 / j23> 

" \ e I: 

, + ‘ 323 n ] 


and the total moment about the parallelepiped's mass center is given by 

- - A 1 31 l 32 l 33 p« [g 312 cos n 312 ^ 312 ^ + g 313 cos 0 313 {^ 13 } J 
+ G 323 cos o 323 ^ 2 ^ ] {«« + {fi 3 ^ {F3 f } d (10 > 

where G^ a ^, H 3a ^ pertain to those sides of area and 


cos n Jae -= { eJ “ e }j T { 6f }j (11) 

Once again, is the vector position of the geometric center of the 
parallelepiped as measured from its mass center. 
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3. 


In the previous section, force and moment expressions were given for 
the four basic allowable shapes. These expressions contained quantities 
which depend on the nature of the flow, quantities which will now be 
defined. 

3.1 Solar Radiation Pressure 

If solar radiation pressure forces and moments are required, then: 

f S Q 

I) P ** P , the solar radiation pressure constant (P s varies as the 

inverse of the square of the distance from the vehicle to the sun, 

-7 2 

and is approximately equal tol.005 x 10 dyne/cm in the vicinity 
of the earth's orbit). 

II) ' a unit vector directed from the sun toward the 
earth. In particular, 

■ W H e < i3 > 


with 


n ■ 


[-sin 0 cos 0.13tt 
sin 0 s sin 0.13tt 


-cos 0‘ 


where 


365.24 


D [radians] 


g 

and D is equal to the number of days after the autumnal equinox 
(it is assumed that 0 s is a constant for the period being considered 
and D S is a mean value). In Equation (13), { 6 } is obtained from 


the auxiliary calculations. 
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Ill) 


(16) 


G J = l-v J 
f = 2v j 

where v~* is the solar reflectivity coefficient for the appropriate 
surface. 

For the sphere and flat plate, only one value for need be 
specified. 

For the cylinder, two distinct values for may be specified, 
one for the side and one for the ends of the cylinder. 

For the parallelepiped, three distinct values for may be 
specified, one for each pair of parallel sides. 

3 . 2 Ae rodynamic Pressure 

If aerodynamic pressure forces and moments are required, then: 

f cl 

I) P = P , twice the dynamic pressure. In particular, 

P a = p V 2 (18) 

2 

where V is the square of the magnitude of the spacecraft inertial 
velocity. 


, 2 . M 1 {r 1 } 


2 —9 

in units of [L 'T ] 


Here, {& r } e is obtained from the Orbital Subroutine and p is the 
atmospheric density.. 

II) - {^ S }j > a Un ifc vector along the negatiye of the spacecraft's 

inertial velocity vector. In particular, 

{ 4 ‘}j “ - Ml Me (20) 

where lb 17 } is obtained from the Orbital Subroutine. 
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Ill) G^ f = a j 

H jf = 2(l-o^) 

where is the input aerodynamic reflectivity coefficient for the 
appropriate surface. 

For the sphere and flat plate, only one value for need be 
specified. 

For the cylinder, two distinct values for a 2 may be specified, one 
for the side and one for the ends of the cylinder. 

For the parallelepiped, three distinct values for may be specified, 
one for each pair of parallel sides. 


3. 3 Gravitational Force and Torque 

I) The gravitational force acting on the rigid Body j is given by 

{f 3 ®}. =-^ 3 [{WM r - { A re }{a r } )] 

3 <'Ti r 'T r 6 e 


( 21 ) 


where y is the input gravitational constant for the orbited body as used in the 
Orbital Subroutine, is obtained from the Orbital Subroutine, and m"^ is the 

mass of Body j . 

II) The gravitational torque on Body j with respect to its center of mass is 


{t JG } , {a 3 ®} {a r } { A j 4{a r } 

J (R r ) 3 6 j 


( 22 ) 


where 


IW. 


are the components in the Body j axes of the centroidal inertia 


tensor of Body j . 
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3.4 Magnetic Torque 


I) The magnetic torque on Body j resulting from its interaction with the 
geomagnetic field is given by the expression 



( 23 ) 


where 



is the magnetic field strength of the earth at Body j and 



is the magnetic moment of Body j. The geomagnetic field strength at 

(*) 

Body j can be approximated by the expression 



3M e 

(R r ) 3 


sinX^ sin\2 COSX 2 

sin^X2 - 1/3 
cosX^ sinX 2 cosX ? 


I 


(24) 


Where M e is the magnitude, of the earth's magnetic dipole moment. Typically, 

o c ■ O 

M e = 8.06 x 10^ oersted-an for the earth. Expressing equation (24) in 
Body j coordinates gives 

II) In the above expressions 

X 2 = sin' 1 {d 2 F {g MG } {d G s } { A er | {d 3 } 


(24.1) 


X^ = sin 


- 1 ' 


KP M 


kl 

W 


M 

H 

CN| 

Q 

{ci 


H 



(. 26 ) 


(*) 


Equation (24) is based on a tilted dipole model of the geomagnetic field. 
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4. SPECIFICATION OF INPUT AND OUTPUT QUANTITIES 

The required input and output quantities will first be specified for 
the case of solar radiation pressure disturbances. Conventionally, (IC) 
denotes an input constant and [1] denotes a non-dimensional quantity. 

4.1 Solar Radiation Pressure Disturbances 


General Input Quantities 


P s 

9 

solar radiation pressure constant (IC) 

[ML 1 T 

D S 

9 

mean number of days after autumnal equinox (IC) 

[Days] 


} • 

transformation matrix 

[1] 

4.1.1 

Sphere 



R js 

> 

radius of the sphere (IC) 

[L] 


9 

position vector of the sphere's geometric 
center as measured from its mass center (IC) 

[L] 


9 

solar reflectivity coefficient (IC) 

[1] 

4.1.2 

Flat Plate 



aJ 

> 

area of the flat plate (IC) 

[L 2 ] 


9 

i 

unit vector normal to flat plate (IC) 

t 

[1] 

M 

j 

position vector of the plate's geometric center 
as measured from its mass center (IC) 

[L] 

j 

9 

solar reflectivity coefficient (IC) 

[l'J 

4.1.3 

Cylinder 

i 


R jc 

9 

radius of the cylinder (IC) 

[L] 

h Jc 

9 

height of the cylinder (IC) 

[L] 

pn 

9 

unit vector along axis of the cylinder (IC) 

[1] 


j 
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position vector of the cylinder's geometric 
center as measured from its mass center (IC) 



solar reflectivity coefficients for sides and 
ends of the cylinder (IC) 


4.1.4 


Parallelepiped 


[L] 

in 


L^W 3 


{e 31 % {e J13 } r , 

W, . 


lengths of sides (IC) 

unit vectors normal to sides (IC) 


position vector of the parallelepiped's 
geometric center as measured from its mass 
center (IC) 


112 il3 123 

r , v J , v J , solar reflectivity coefficients (IC) 


M. 

H 


General Output Quantities 
, solar force on Body j 

, solar moment about Body j mass center 


[L] 

[1] 

tL] 

tl] 


[MLT -2 ] 

[ml 2 t -2 ] 


4.2 Aerodynamic Pressure Disturbances 


H 

M 

H 


General Input Quantities 


velocity vector of the orbital reference frame 

[LT _i ] 

[ 1 ] 

[ 1 ] 


transformation matrix 
unit vector parallel to R 


Sphere 


R JS 


radius of the sphere (IC) 



position vector of the sphere's geometric 
center as measured from its mass center (IC) 


0 


j 


aerodynamic reflection coefficient (IC) 


[L] 
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[L] 

[13 



4.2.2 


Flat Plate 


aJ 

> 

area of the flat plate (IC) 

[L 2 ] 

M, 

9 

unit vector normal to flat plate (IC) 

[1] 


9 

position vector of the plate’s geometric 
center as measured from its mass centeir (IC) 

[L] 

o j 

9 

aerodynamic reflection coefficient (IC) 

[1] 

4.2.3 

Cylinder 



R jc 

9 

radius of the cylinder (IC) 

[L] 

h jc 

9 

height of the cylinder (IC) 

[L] 

M 

i 

9 

unit vector along axis of the cylinder (IC) 

[1] 

H 

9 

position vector of the cylinder’s geometric 
center as measured from its mass center (IC) 

[L] 

is ie 

a J , a J 

9 

aerodynamic reflection coefficients for sides 
and ends of cylinder (IC) 

[1] 

4.2.4 

Parallelepiped 


i 11 . l^ 2 . 

3 , 

lengths of sides (IC) 

[L] 

|e^ 2 |. . ,je^ 2 ^| , unit vectors normal to sides • 

^ ^ i 

[1] 



position vector of the parallelepiped's 



geometric center as measured from its 
mass center (IC) 

[L] 

c^V 13 , 

j 23 
.0 > 

aerodynamic reflection coefficients (IC) 

[1] 


General Output Quantities 


l F j* 

1 : 

aerodynamic force on Body j 

[MLT -2 ] 

{ T jA 

aerodynamic moment about Body j mass center 

2 -2 

[ML T ] 
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4. 3 Gravitational Force and Torque 


General 

Input 

Data 


Y 

9 


gravitational constant for the orbited 
body (IC) 

[l 3 t“ 2 ] 

m^ 

9 


mass of Body j (IC) 

CM] 

R r 

9 


magnitude of orbital radius vector 

[L] 

Ml 

9 


transformation matrix 

tl] 

Mr 

9 


vector from orbital reference to Body 
mass center 

j 

[L] 

M} 

9 


transformation matrix 

[1] 

{•1. 

9 


— ■£* 

unit vector parallel to R 

[1] 

ki 

9 


transformation matrix 

[1] 

M 

9 


Body j inertia matrix (IC) 

[ML 2 ] 

General 

Output Quantities 


W. 

W, 

9 


gravitational force on Body j 

[MLT -2 ] 

9 


gravitational torque on Body j 

[ml 2 t -2 ] 

Magnetic 

Torque 

• 


General 

Input 

Data 


Mi 

) 


magnetic moment of Body j (IC) 

[ML 2 T 2 (oersted) 1 ] 

M e 

9 


magnitude of earth's magnetic 
field (IC) 

[L J (oersted)] 

U 3e } 

• 


transformation matrix 

[1] 


[HL 2 T _2 (oersted) _:L ] is equivalent to [L 3 - (oersted) ] . The latter 
form is included for the user's convenience, since the torque then 
has the proper dimensions of 

. jM, _ [ML 2 T~ 2 (oersted)" 1 ] [L 3 - (oersted)] » ML 2 T 2 

J = O 
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A re j , transformation matrix [1] 

t G , mean time in seconds after noon (GMT)(IC) [sec] 

D S , mean time in days after autumnal equinox (IC) [days] 


General Output Quantities 


{A 

Hr 


magnetic torque on Body j 

earth's magnetic field strength at Body; 


9-2 _ 1 

[ML T (oersted) ] 

[oersted] 
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6. 3 Flexible Body Environmental Disturbance Equations 


The following equations constitute the generalized forces acting 
on the flexible bodies of the spacecraft model due to solar radiation pressure, 
aerodynamic pressure and gravity gradient effects. Derivations are con- 
tained in Appendix B. 

Note here that those environmental disturbances of concern in the 
analysis of the Skylab vehicle have all been included (solar, aerodynamic, 
gravity) while magnetic disturbances are neglected since the design of the 
various flexible arrays (negligible dipole moment requirement) precludes the 
problem of magnetic interaction. 


R^ S 

a 

. P S A 3 A 3b<T) { V 3 

b 3 o (B 

j(T) 

6 

A 3 / S e 
6Y Y 

) 



+ (1 - p3 > A is s e } 1 

B 3 <T) 

a 

A je 

ap 

s e | 

p 1 




■ pS a3 *L{ vJ 4 


A je 

6Y 

s e 

Y 




+ a - pJ > a sb s e } 1 

b KT) 

A je 

ap 

s ;i 




■ pS aJ + “is 

){ V 3 

D j 
B 6 


A je 

<$Y 

b r 

Y 


+ (1 - ^ A j« S 5 } 

|b 3(i) 
1 a 

A je 

ap 

B pi 



R jA 

a 

- -pV 2 A 3 A^ (Ti { (1 

-ah 

b b 

<4 m 

A je 

<$Y 



+ A SS b 6 } l B o <T) 

A 36 

ap 
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<4* 

■ - (,v2 aJ *kB { (1 

-o 3 > 

B & 

<4 m 

A je 

Sr 



+ * 3 41 

a 3 ‘ 

ap 
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( 6 - 13 ) 


( 6 - 14 ) 


( 6-15 


( 6 - 16 ) 


06 - 17 ) 
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- - pv 2 aJ <;° b + H 33 ) { (1 ' <^> B 3 (B j s (I> A 2 * b') 


A ^ 6 S 6 

A g 6 b 6 


} i*J w AJ_ b 


i( T ) A J e b r i 

<jp p 1 


JG 


ynr 


(R r ) 3 


a 3 


JG = _ ym 
K (R 


[4 - 3 4; *; < 4 <T) *n •? 

4 3 >] 


+ A^< T) h| 3 - 3 A- (a 


'oB 3 


^Je(T) 

Pp 


— r 

r > 3 L 


♦j A ift *{ + °’ j + 

ka c*3 3 


4 


"it 4 - 3 4 a 4 
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2 HE 
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ap 


A j£ 
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A 

A 3e<T) 

I jf 

A je 

a6 

36 

6e 

ecf 

pp 

s j6 

a3 

A jr 

36 

n 6 ‘ 3 

s j6 

a3 

A je 

A 36 


r a ](T) A re & r ) 


3 6 (n p 


(R )' 
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(R r ) 3 L 

+ 3A J« > A je(T) (bJ (T) + S^ 1CT) + S^ 3(T) ) qjl 

+ 3 a 3 36 6 s ^ ek ek eR k J 


Pp 


(6-18) 


(6-19) 


(6-20) 


( 6 - 21 ) 


where 


n a " 

R 1 + 
a 

d.r (T) 
a3 

n ij 

3 

(6-22) 

D k " 

7 f .< 

dm 3 

k 

(6-23) 

b L = 

r (T) 
a 

e 

je (T ) 

e3 

B j 

B k3ct 

(6-24) 

Bp » = 

k6s 

hf 3 

m b j 

r 3 ({i? dm 3 

k 

(6-25) 

i ■ 

r (T) 
a 

e 

A 3e(T) 

ea 

C j A je a r 

kJIoB 36 6 

(6-26) 

r 1 - 

k£aS' 

m J B- 

P 

| \ 

^ dm j 

(6-27) 
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(6-28) 


S 

S 

s 


jll 

ka 

J12 

kJt$ 

jl3 

la 


a 


a 


q 


r(T) 

Je(T) 

B;| 

ka0 

S 

% 

r(T) 

A je(X) 


e 

ea 

k£.af5 


j c j 3.2 

k k£a 


(6-29) 

(6-30) 
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VII. Orbital Subroutine 


The dynamic equations for a multibodied flexible spacecraft presented 
in the previous sections of this report are written with respect to a 
reference axis frame. One important choice for specification of this reference 
axis frame is that it lie on a user-defined Kepler orbit. The orbit routine, 
as specified in this section, calculates this reference orbit. The orbit 
equations contained herein were specified concurrent to the rigid-body 
environmental disturbance equations and thus are likewise included here 
in their original form in order to retain correspondence with the computer 
programming documentation. As in the case of Section 6.2, matrix notation 
is used herein instead of the index notation employed throughout the body 
of this report. (The orbit subroutine commences here with Section 2.0 
following. ) 
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2.0 COORDINATE SYSTEMS AND NOTATION 


In describing the motion of the reference axis frame, two right-handed, 
orthogonal coordinate systems are of basic importance. These coordinate systems 
are defined and described below. The inertial reference frame is defined with 
respect to the earth as the orbited body. In case orbital motion about another 
celestial body is desired, an appropriate inertial reference frame must be defined 


2. 1 Inertially Fixed Coordinate System 

6 6 6 

The inertially fixed coordinate system with axes (x^, x x^) has its 
origin at the mass center of the earth and is defined by unit vectors e^, 
e^» directed along the appropriate axes in a positive sense. The unit 
vector e!j normal to the equatorial plane and positive northward, while 
e® is directed along the autumnal equinox. The coordinate matrix for this 
system is given by 



( 2 . 1 ) 


Since the coordinate system is inertially fixed, its rotational velocity vector 
is given by 



( 2 - 2 ) 


2.2 Orbital Reference Coordinate System 

X X T 

The orbital reference coordinate system with coordinate axes (x^, X 2 » x^) 
is defined by unit vectors e_^, e^, e_^ such that e_^ points toward the earth s 

mass center, ef lies in the orbit plane and forms an acute angle with the 

*** ^ ^ 
tangential velocity vector while e^ is normal to the orbit plane. The coordinate 

matrix for this system is given by 
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and its angular velocity vector is denoted by 


(2-3) 


T? 4 o 


(2-4) 


The vector distance from the center of mass of the earth to the origin of 
/ r r r N , :rr 

(x^ x 2 , x 3 ) is R . • 

2. 3 Notation 

The following common notational rules will be used in the sequel: 


(i) The vector f is represented as a column matrix by | f where j 
denotes the coordinate system in which the components of f are 
given; i.e., 



j equal to r implies components in the orbital reference 
• r r r 

coordinate system (x^, x 2 » x.j), while j = e implies components 
in the inertial reference system (x®, x 2 , x^)* 

(ii) A A above a symbol denotes an outer product matrix; i.e.. 
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(iii) A superscript T on a matrix denotes its transpose. 


— — IT * *T 

The following sections are devoted to the determination of R , U) and 
their required time derivatives. 


3 *° ORBIT EQUATIONS 

j* 

Let us now define the orbit equations for the reference frame 
The position of the reference frame is determined by the equation for an inverse- 
square central force 


where 


d 2 R 1 


dt 


F 

M 


F - 
g 


CR r ) 3 


9 


with: 


M «=» total spacecraft mass [M] 

R r = magnitude of "r 1 [L] 

3 -2 

Y = gravitational constant for the orbited body (GM^) [L T ] 
where G is the universal gravitational constant 
and M, is the mass of the orbited body 


(For the earth, y = 3.98604 * 10 5 km 3 /sec 2 - 1.40766 x IQ 16 ft 3 /sec 2 ) 
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Thus, the position of the reference frame is determined by 


d 2 ^ 




(R r ) : 


(3.1) 


Equation (3.1) represents the motion of a point in an elliptical orbit about 
a spherical planet. Figure 3.1. 



Figure 3.1 Orbit Parameters 
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Letting 


P = orbital period [T] 

e = eccentricity of orbit [non-dimensional] 
the semi-major axis of the orbit is given by 

(3-2) 


(3-3) 


fr) 


2/3 


, W 


where the "mean motion" to is defined by 
o 


2 IT r 1 


to 


, (T A ] 


with to Q being the equivalent circular orbital frequency in radians per 
unit of time. 


The true anomaly t|t is determined by 


= 2 tan 



1+ e 
1- e 



(3-4) 


where the eccentric anomaly E must be found as a solution of Kepler's 
equation 


E - s sin E = to (t - t ) . 

o p 


(3-5) 
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Here, t is the time at perigee and the quantity w (t - t ) is the mean 

p op 

anomaly. 

The radial distance to the reference frame, R r , is given by 


? <A- . ‘3 _. [L] 

1 + ecosy 


In addition. 



V 


ayd-e ) 

R r2 



] 


By differentiating the above relationships we find that 


(3-6) 


(3-7) 


R 




rU 


(1-E 2 ) ^ 3±n ’^ ^ LT ^ 


(3-8) 


R r = 


eycos \p 
„r 2 


[LT -2 ] 


(3-9) 


2 eysinijj 
R r 3 


[t" 2 ] 


(3-10) 


4.0 MOTION OF ORBITAL REFERENCE SYSTEM 

6 6 6 

The inertial reference system (x^, X 2 ,x^), the orbital reference system 
(x£, xljjX^), the transformation between them { A r e } and the angular velocity 
vector 03 of the orbital reference frame are independent of the relative motic 
of the bodies of the spacecraft and can be written as functions of the orbit 
parameters and time.. We now determine these functional relationships. 

4.1 The Inertial Reference Axes 

The inertial reference axes have their origin at the center of mass of 
the earth. 
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0 

The axis is normal to the equatorial plane and is positive in the north- 
0 

ward direction; the x., axis lies along the autumnal equinox. 



Figure 4.1. Inertial Reference System 


The coordinate matrix of the inertial axes is given by (2-1) as 



4.2 The Orbital Reference Axes 


p 

A rotation of the inertial axes through an angle S about the e~ axis 

/ b! 8 ^ 

produces the \ efysystem in which lies along the line of nodes (intersection 
of the orbit plane and the equatorial plane, Figure 4.2) with e^ directed toward 
the ascending node. 
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Figure 4.2. Line of Nodes 


With 


e P 

1 

1 

let S A^ e f be the transformation matrix which 


l “ / 

—2 


l > 


transforms the | e 6 ^ system into the |e^| system such that 

K} ■ K1U1 

where 



(4-1) 
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and 


-sin/3 0 -cos/3 

00 0 (4-2) 

cos/3 0 -sin/3 - 

Rotation of the ) e^ f system, Figure 4.2, about e& through an angle 4, the orbit 

l - J , t ) _3 t 

inclination angle , produces the < er > system. Figure 4.3, where ej lies in the 

orbit plane. 




Figure 4.3. Orbit Inclination 
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Since £ is a constant, £ 


0 and 



(4-4) 


Rotation of the je^ system. Figure 4.3, about e| through an angle a produces the 
system. Figure 4.4, where e^| passes through the orbital position of the 
reference frame. 



Figure 4.4 Orbit Angles 
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With 




let 



be the transformation matrix which transforms 


the 




system into the j e^^system such that 


K}- 

where 

cosa 0 -sina 

0 10 

sina 0 cosa _ 

and 

-sina 0 -cosa 

0 0 0 
cosa 0 -sina. 




(4-5) 


(4-6) 


The orbital reference system (x*, x^, x^) with unit vectors e^, e!j» 

£3 has its origin at the orbital reference point, with eT pointing toward 
nadir, e^ in the orbit plane forming an acute angle with the velocity vector 
and normal to the orbit plane, Figure 4.5. 
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Figure 4.5. Orbital Reference System 


Thus, 



Letting 




and letting 



be the transformation matrix 


which transforms 


the 


{^system into the|e^ r 


system. 


we have 




Combining the previous equations , 


or, letting 


where 


and 


Similarly, 


or letting. 


KH A i{ A “ 5 }{4 



-sindf 0 


-cosa 


( 4 - 9 ) 


{A r ^}= a 


-sina 

0 

-cosa 


0 

0 

0 


-cosa 

0 

sina 


( 4 - 10 ) 


{e 1 }- {a^WIM 
H - {/I . 

{ 4 = { a "}{4 <4 ' u> 


77 



where 


w = 


(coso;cos|cos/3-sina;sinj3) (cosasinl) (-cosacos^sin/3-sin«cos/3) 


sin|cos/3 


) ( -cos| ) 


-sin| sin/3 ) 


(-sinarcos/3cos|-cosa:sin/3) (-sinasin|) (sinacos£sin/3-cosacos/3 ) 


Thus, {A re | is a function of the three orbital parameters /3 , land a. 

The angle a has been defined as the orbital position of the reference frame with 
respect to the ascending node while*/' has previously been defined as the true 
anomaly or as the orbital reference position with respect to perigee. Therefore, 
if is the orbit angle between the ascending node and perigee, then 


a - a + ip 

P 


and from Equation (2-7), (2-10): 


&= $ = 


aT (1 -e 
_ 2 


, [t" 1 ] 


(4-13) 


a= ip - 


-Ze.y s in tp 
_ 3 


[T" 2 ] 


(4-14) 


78 



The orbital inclination, ^ , is an input constant. The nodal regression 
rate (3 due to the orbited body's oblateness is given by 



0/1 2.2 
2(1 - e ) 



(0 

o 


cos 5 [T 


(4-15) 


where R^ is the radius of the orbited body and is the 

second harmonic of the orbited body's gravitational potential. For the earth, 
\ 13 6.378 x 10 3 km = 2.093 x 10 7 ft 

J 2 = 1082.28 + 0.3 x 10“ 6 

(The above expression for 0 as well as the value for of the earth are found 

in the Reference below (*) . If 0 is the value of 0 at time of the first 

o 

nodal crossing, t , then the parameter 0 is given by 


0 « fi Q + 4(t-t n ) . (4-16) 

Since the matrix is zero (Equation 4-4), differentiating the relation 

for |A re j gives 

{A re } - {A r *}{A W } {A 06 } + {a^}{a^}{a^} . (4-17) 

But, it can be shown that 

<*-»> 

Rewriting Equation (4-17), we find that 


M - H Kf { {^1 {a p ? kfM 

so that 

-{**«} {a rf } T - {*«}{/»} {/*} trY {*»} W ■ (4-19) 


(*)Escobal, P. R. , Methods of Orbit Determination, John Wiley and Sons, 
Inc., New York, New York, 1965. 
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Expanding the above expression 


where 






sin^cosQ'l 
-cos? I 
1 -sin?sinoi 



► 


Therefore, 



j3sin?cos a 

- 0 cos£-a 

- 0 sin?sintt 


( 4 . 20 ) 


and, since ? is constant. 



• • 

-0a sin ? sine* 
-0a sin?cos a ' 


( 4 - 21 ) 


*» 
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Finally, letting 



( 4 - 22 ) 

( 4 - 23 ) 

( 4 - 24 ) 

( 4 - 25 ) 

( 4 - 25 . 1 ) 

( 4 - 26 ) 

( 4 - 27 ) 
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5.0 CONCLUSIONS 


We have here documented the determination of those orbital parameters 
necessary to specify the reference axis frame in the given Kepler orbit. 

In particular, the requisite input parameters are simply: 



= orbit angle from line of nodes at time zero, [radian] 

= orbit angle between ascending node and perigee, [radian] 

=* angle between autumnal equinox and the ascending line of nodes 
measured in the equatorial plane at time t R . [radian] 

= orbital period 1^] 

= orbital eccentricity. [-} 

= orbit inclination with respect to equatorial plane. [radians] 

= radius of the orbited body [L] 

3 -2 

= gravitational constant for the orbited body [1* T ] 

= second harmonic of the orbited body's gravitational potential [-] 
= time of first nodal crossing (-P^t^iO) [T] 


The output parameters are: 


R 

V 


Mr 

fir 

w 

He 

He 

ft 

He 

M. 


magnitude of orbital radius vector R [L] 

-r -I 

magnitude of orbital velocity vector R [LT J 


.-li 


inertial angular velocity of orbital reference frame [T J 


„- 2 1 


inertial angular acceleration of orbital reference frame [T ] 

transformation matrix from inertial to orbital reference axes [ 
column vector representation of R . [L] 

unit vector along R . [-] 

* 'y -.1 

column vector representation of R . [LT ] 
unit vector along R . [-] 

.Hr -2 

column vector representation of R . [LT ] 
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6.0 EXISTING COMPUTER PROGRAMS 

Since a digital orbital routine identical to the one presented in this 
report has been successfully coded in the existing TRW Generalized Spacecraft 
Simulation (GSS) program, this routine is included here for 

reference purposes. (Note change in expression for 0 and inclusion of (59)- (65).) 

6.1 Subroutine SET0RB 

This subroutine initializes the orbit parameters and computes the con- 
stant portion of the orbital computations. The initial eccentric anomaly (EK) 
is computed and used as the starting value for the first iteration of Kepler’s 
equation. TPER, the time at which the spacecraft last crossed perigee is com- 
puted. The semi-major axis of the orbit (AXIS) and a related constant (AESQ) 
are computed from which the nodal regression rate (BD0T) is computed. 


SUBR0UTINE SET0RB 


a 

<P 


\ 



sin (ip) 
l.+cos (ip) 


* TPER - (E k esin(E k ))P 


2ir 


* a 
AESQ 


t. 


(p/ 7/2,) 2/3 
a(l.- e 2 ) 

3 -. 


2. (1-e 2 ) 2 


(-H 


2tt cos (Q 
P 


83 



6.2 Subroutine 0RBIT 


This routine computes the orbit angle ALPH (the true anomaly measured 
at the center of the earth between perigee and the origin of the reference 
frame) at time T by solving Kepler's equation using Newton's method of 
iteration. Each iteration is begun by assuming the previous value of ALPH 
as a starting value. A maximum of ten iterations is allowed, the iterations 
ending when the difference between two successive values is less than 10 
If this condition is not met within ten iterations, the initial value is 
changed, a message is printed out and the program begins the next case. 

All orbit related parameters are computed within this routine, including 
the orbit radius (RR),tlie radial velocity (RRD0T) and acceleration (RRDDT) , 
the orbital angular velocity (PSD0T) and acceleration (PSDDT) , the angular 
velocity of the orbital reference frame relative to inertial axes (0MRRC) 
and its derivative (0MRRD) , the gravitational acceleration at the reference 
frame origin (C0NG) , and the angle between the autumnal equinox and the line 
of nodes (BETA) . 


SUBR0UTINE 0RBIT 

* ( 1 ) 

(5) 

(7) 


( 8 ) 

( 22 ) 

where KP = Number of Restarts 


ECOUT = 2ir (T — TPER)/P 
SL0PE = 1. - ecos (E^) 

STEP = \ ~ esin (E^ - ECOUT 
SL0PE 

Ejj = \ “ STEP 

E fc = Ej^ (1. - 0.1 KP) 
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(36) * 

(44) ECPS 

(45) R r 

* (46) C0NG 


2. T 


-1 


an 


(h) 

* l.-e an \ 2 . / 


1. + e cos (ip) 

AESQ/ECPS 

2 

y/R r 


* (47) 

(4 8) 

(49) 

(50) 
(53) 

(57) 


(58) 


R r 




4> 

P 


R 


M 


e sin(t}>) ^/y/AESQ 

C0NG e cos (0) 

'2 


YVAESQ /R r " 


-2 . V e 


V “Ct-O 


sin(i/))/R r 


0 s in (4) cos(d;) 

~0 cos (4) -‘Z 1 
-0 sin (4) sin(a ) 

-00 sin(4) sin(a) 

-ip 

-0 ip sin(4)cos (Q!) 


(59) 


M - 


(cosacos4cos0-sinasin0) (coso;sin4) 

( sin4cos0 ) ( -cos4 ) 

(-sina!cos/3cos4-coso!sin0) (-sinasin4) 


(-cosacos4sin0-sinu?cos0) 
( -sin4 sin/3 ) 
(sinacos4sin0-cosacos0 ) 
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* ( 60 ) 


M- ° 

( 1 1 

* (61) {r 1 } = -R r {a"} T {s r } 

e 

* < 62 > M e ■ 

* (63) {tv 

(63.1) V 2 - {R r }^ {»'}„ 

* < 64 > Me - Me /V 

* (65) {ii r } e - - {a 16 } (it r {s r }+ 2 R r M r |s r } 

♦ R r {;*} r H * * r M c M r M) 
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kp - kp*i 





ra 

213J 

< 

KP: 


ER0RB 

■ 1# 


Slope of Kepler's Equation. 


Zero slope Is error condition. 


Compute one Iteration of Kepler's 
equation using Nevton's method. 


Equation Is solved when Iteration 
Is less than 10*5. 


Iterate 10 times. 


If do solution, restart vlth 
different initial value. 


Restart k times maul as s. 


If no solution, exit vlth 
error flsg set. 


( RETURU 3 
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Compute orbit angle 


Compute orbit radius 


Gravitational constant at orbit 
radios and radial velocity. 


Radial acceleration and angular 
velocity. 


Angular acceleration 


Angular rate of orbital axes in 
orbital coordinates. 


Angular acceleration of orbital axes 
in orbital coordinates. 














VIII. Control Subroutine 


The UFSS Program contains a generalized control subroutine capability 
which allows the user to quickly and efficiently synthesize almost any con- 
ventional continuous attitude control system. The control system simulation 
capability is fully described in Reference [4]; this preent section describes 
the interfaces between the dynamics and the controls subroutines. 

Figure 8.1 presents a data flow diagram showing the various features 
of the Dynamics /Control Generalized Interface. The Sensor Interface generates 
the transformation matrix and angular velocity components 

Ctp 

necessary to simulate the nth attitude and/or rate sensor mounted at node 
N (or at position r*^ ) on a given Body j of the spacecraft model; these 
quantities are then transmitted to the Control Routine for use in performing 
the attitude and rate error calculations. 

The Generalized Force Interface accepts as inputs the forces and 
torques produced by the Control Routine and transforms them into the gen- 
eralized forces required by the Dynamics Routine. 

If the control and steering laws (or any relationships existing in 
the control system) involve the solution of a differential equation, then 
the appropriate dependent variables are referred to as components of a control 
state vector. These control state vector equations must be added to the 
dynamic state vector equations and integrated simultaneously; this combining 
of the dynamic and control state vector equations is effected in the State 
Vector Augmentation Interface. 

8.1 Sensor Interface 

Assume that the nth control element (sensor or controller) for Body j 

is located at node N . In order to locate the elements sensitive axis, it 

is necessary to determine the instantaneous orientation of an axis frame 

IN 

located at node N and fixed in Body j. Specifically, if r is an axis 

j ^ 

frame based at node N and parallel to the er frame for zero flexible 
motion of Body j , then it is necessary to determine the time-dependent trans- 
formation matrix A^ relating e^ N to e| for Body j in its perturbed 

ap -a -p 

position: 
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where in general is a function of the q^ (£=1,2 n^). 

When Body j undergoes flexural motion, the axis frame rotates 

1 -a 

from the e^ orientation to its final orientation as shown in Figure 8.2. 




To a linear approximation, this rotation of the e- 1 frame is expressed 

-a , 

by specifying three small rotations about an axis set e J parallel to the 

e^ set but located at node N . Calling these rotations , it follows 

u oi 

that 


4 > 


JN 

a 


ak 


e-% q 


( 8 - 1 ) 
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JN 


e-^' is the 3 x n. matrix whose kth column contains the three modal 
ak J 


where 

rotation components at node 

JN 


Since the ijr 


N of Body j for its kth mode, 
are small rotations, it follows that 


e** 

-a 


= A' 


jNj 

ctS 


4 


where 



1 

JN 

*3 

JN 

”*2 

A jNj - 
aS 

,jN 

-*3 

1 

JN 

n 


JN 
- 2 

JN 

“*1 

i 

or more simply: 

a$ 

= 5 aS 

- « 


( 8 - 2 ) 


where 



is the unit 'tensor or identity matrix. 


Now let us locate the nth control 

it e^ cn , at node N . The frame e^ cn 
. . -oi -a 

e^ N by a constant input transformation 


element reference frame, call 


is related to the body-fixed frame 
. jcnjN . 
aB 


where 


e 3 cn = A JcnjN g jN 
-a aS -8 


A 3 J nJ " - (G 3< =") . (G 3cn ) 


08 


ay 


ya 


as 


( 8 - 3 ) 


* i on 

The (G" 1011 ) „ are formed from three user supplied Euler rotations 8“ as 
v p aS P 

follows: 
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1 

0 



0 

(G jcn ) 

^1 

0 

cos 

e jc» 

sin e| cn 


0 

-sin 

„jcn 

6 i 

_ jen 
cos 0^ 


cos 

Jen 

0 2 


0 


/Jen. 
(G 2 } a8 

0 



1 

0 


sin 

Jen 

9 2 


0 

Jen 
cos 0^ 



.i en 
cos 0^ 

sin 0f n 

r 1 

o 

, Jen 
(G 3 } a8 

. Jen 
-sin 0^ 

cos 

0 


o 

0 

1 

J 


(8-4) 


(8-5) 


( 8 - 6 ) 


i cn 

Finally then, the nth control element reference frame e“* is related 

to the Body j reference axis frame e^ through the time-dependent trans- 
, . j cn j 

formation matrix A^ 


jcn 

e J 

-a 


= A jcnj 

a 6 -8 


where 


A jcnj = A jcnjN A JNj > 

a$ ay y8 


that 


J ^ j 

It should be noted that if Body j is rigid, then A 

A J cn j _ A jcnjN , a constant matrix. 
aB a8 



(8-7) 


) 
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In addition to calculating the orientation matrix 

— -I QJJ 

Interface calculates the inertial angular velocity or 

of the e^ Cn frame, where 
-a 


A jcnj 
= ur 


» 

cn 


the Sensor 

icn 

e 


-jcn -j 
or = or 


+ 


with 0*^^ being the angular velocity of the e^ Cn ftame relative to the 
e] frame. In particular, 

-p 

a -a 


where 


;,jN _ 


a 


JN n* 

6 ak q k 


Thus, 


4 C " - A irs + H" 1 


8.2 Generalized Force Interface 

If the nth control element (located at node N) for Body j is a 
thruster, then the output of the Control Routine is the instantaneous value 
of the thrust vector components in the nth control element reference frame, 

pint = F jnt e jcn _ 
a -a 

If the nth control element is a torquer (CMG) , then the output of 
the Control Routine is the instantaneous value of the torque vector com- 
ponents in the nth control element reference frame, 

ifjng = T jng e jcn 
ot -a 

If the nth control element is a gimbal hinge torquer, the output 
of the Control Routine is the instantaneous value of the motor torque 
components about the hinge axes, 


( 8 - 8 ) 


(8-9) 


( 8 - 10 ) 


( 8 - 11 ) 
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If Body 1 is a rigid body , then the output of the Generalized Force 
Interface is the total force and moment about the mass center 



Lj ctg g 

n 




A jcnj(T) 

ag 


T j n S . Ijcn jcnj(T) 
T g + r ctg A gy 



where r^ cn = r^ Cn e^ is the input position vector of the nth control 
a -a 

element on Body j with respect to its mass center. 


If Body _j is a flexible body , then the Generalized Force Interface 
must calculate the required generalized forces. Let us first assume that 
the nth control element is a thruster at node N . 

The virtual work done by the thruster force df^ tn acting on the 
arbitrary mass element dm^ of Body j due to the virtual displacement 

6R 1 + 6^ + <$6^ x (r^ + u^ ) 


6W^ nT = J df^ tn . JjSR 1 + <Su^ + fig'’ x (r^ + u^)J 

B^ 

where 


6R 1 

= 5R 1 

r 

e 


a 

-a 

<5u^ 

= 6q j 



= 60 ^ 

e^§ 


a -a 


But, for a thruster mounted at node N , 

d f j tn = F jtn 6 - ? jN ) d V j 


( 8 - 12 ) 

( 8 - 13 ) 


( 8 - 14 ) 
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where 


pjtn 



dV 


j 


= thrust vector from the nth thruster on Body j located 

i cn 

at node N (in e J frame) 

-a 

= position vector of node N relative to the Body j 
hinge (in e^ frame) 

= arbitrary volume element of Body j 




6 (a - a ) 
o 


Thus 


position vector to an arbitrary point within dV 
relative to the Body j hinge (in e^ frame) 

Dirac delta function having the property that 

f f(a) 6 
J 


j 


(l - a ) dV j = f (a ) 
v o o 


J U > T , fjtn | 5r 1 s r 


+ S< 4 4 + 


60^ e^ 8 
a -a 


,-jN . 

x (r + 


; JN )] 


or 


6W JnT = R^ nT 6R 1 + Qp nT 6q^ 

a a k k 

so that 



R^ T = 


r 

e 


a 


-a 


r - 

F Jnt . 

e^ 

-a 

JN 

^ka 

Q®r - 

jjnt _ 

[4* 

x + ;J")] 


( 8 - 15 ) 


where R^ nT , > Q(H)^ nT are the thruster generalized forces associated 

ot k ot ^ § • 

with the generalized coordinates R 1 , and 0^ respectively. 

Ot K Ot 
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Thus , 


or 


Likewise, 


so that 


Finally, 


Q® J a 


jnT _ 



rJ- t 

a 

- J nt 
' F e 

Jen 
-6 * 
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e 

-a 


R jnT 

a 

= A rj 

ay 

Jcnj (T) 
Y0 

Jnt 

F 0 



- 4 nt 

jen 
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J JN 

-a ^ka 
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Jcnj<T) 

a0 

Jnt 

f s 

- pj nt 
■ F 0 

j° n 
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■ {4 s = 

/ jN 
« ( r Y 

+ u y N) ~y 

- e jg 
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k N 

, JN, 
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) e^ x 
-Y 

4 nt 4“} 

_ r J(T) 
" G ao 

MK 
( ay 

. ~jN 
+ u 

ay 

l » j cn j (T) jnt 

/ y 0 e 


But, transforming from hinge axis to Body j axes, 

-1 




jnT 

0 


so that 


®r - k + 4r 3(I> 4 


jnt 

0 


( 8 - 16 ) 


( 8 - 17 ) 


( 8 - 18 ) 
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Let us now assume that the nth control element is a torquer (CMG) 
at node N . In order to calculate the generalized forces due to an instan- 
taneous torque T^ n produced by this torque generator, one introduces 
two fictitious nodes N + and N located at the terminus of the vectors 
1/2 e and - 1/2 e respectively from node N . At node a force F^ n+ 

is introduced and at node N a force F^ rt is introduced such that: 

(see Figure 8.3). 


1) F jn_ = - F jn+ 

2) T 3(3n =■ 1/2 e x F 3n+ + (- 1/2 6) X F 3 °~ 

- 1/2 6 X <F 3 " + - F 3 "-) 



Figure 8.3. Introduction of Fictitious Nodes 


In this case. 




(F jN+ 

+ f^) . 

6R 1 e r 


a -a 


fJ n+ . 


4 - : 

F jN - . 5q{ 

F jN+ . 

60^ e^ g x 

( ^ K+ 

. -jN+x 

+ u J ) 


a -a 

F jN " . 

60^ e^ g x 
a -a 


-jN-x 
+ u J ) 



-a 
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or 


6W' 


jnG _ 


R i” G SR a + Q k" G 4 «k + «®f° S6 t • 


( 8 - 19 ) 


so that 


R^ nG - 0 (since F jN+ + F jN 3 0) 

a 


(8-20) 


JnG _ 
‘k 

F jN+ . 

J N+ + F JN ' 

ka -a , . 

Jn- 

^ka 

* 0^ 
-a 

jnG _ 
a 

pjN+ . 

e^ g x (^ N+ + 

-a 

— jN+ s 
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+ F jN " . e^ g 
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+ e, 

k 

1 - 
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fJ"-)] 
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H' ■ 

fjng 




,-jN- . -jN-v 

x (r J + u J ) 


gjN 

6 k 


t- \ o] 


Therefore 


JnG = jN(T) .jcnj(T) jnG 
ka a 8 8 


( 8 - 21 ) 


Now let us examine Q(H)^ nG 
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q dP^ nG 
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e jg 
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F jN+ . 

e^ S 
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Q®r 

= 
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e jg . 
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(q i 

5f 

x | e) 

x fJ N+ 



- 

e jg . 
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H" 

x | S) 

x F jN_ 



a 

e Jg . 
-a 


+ 

e jg . 
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Retaining only 
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order term. 




q @J» g . eJs . i J " 8 

_ pj(T) ,jcnj (T) jng 

' G a$ % Y 


so that 


Finally then. 
Body j are: 



the total control 


jcnj (T) 


generalized 


T jng 

e 

forces acting on the flexible 


r^ c - 

a 

E 

n 

r^ t 

a 


<4 C ■ 

■ E 

n 

< r 

* <T> 

®f ■ 

■ E 

n 

(®i nT 

+ ® J a nG > 


( 8 - 22 ) 

( 8 - 23 ) 

( 8 - 24 ) 

( 8 - 25 ) 
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IX- Hass Properties Subroutine 

The UFSS Program requires that the flexible deformation characteristics 
of the terminal bodies in the spacecraft model be expressed in terms of 
three-dimensional orthogonal functions. The standard source of these 
orthogonal functions is a Structural Dynamics Program (SDP) wherein each 
designated flexible body is modeled in the traditional structural dynamics 
sense as a lumped parameter system. The basic requirements placed on the 
structural dynamic model for a given flexible body are as follows: 

1. The structural dynamics model for each simulated flexible body 
must basically consist of point masses only, with no significant 
lumped inertias admissible (no large lumped inertias are admissible 
because of the manner in which certain mass integrals over the 
flexible bodies are presently calculated) • It is felt that this 
restriction is not a serious one particularly in view of the 
significant increase in calculations and cost necessary to allow 
for large lumped inertias in the structural dynamics model . 

2. All three components of the lumped mass at a given node must be 
identical. 

3. If sensors and/or thrusters are mounted on a simulated flexible 
body, then a node point (joint) must be positioned at each sensor 
and/or thruster location and modal rotations must be available 

at these node points. 

The UFSS Program is designed so that choice of a structures program 
is completely arbitrary so long as a basic data set is obtainable. In 
particular, the SDP provides the following data set for each flexible body: 

1. Model Description 

1.1 joint coordinates (r^) 

J a 

1.2 joint mass values 
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2. Orthonormal Mode Specification 


2.1 eigenvalues associated with each mode 

2.2 joint deflections and rotations for each mode (0^ and 0^) 

It is important to note that no development work was intended on the 
SDP. Depending on the particular choice of a structures program, it is 
only necessary to appropriately reformat the output to agree with the format 
of the Structures Tape as given in the UFSSP Users Manual. 

The Mass Properties Routine utilizes the modal and structural dynamics 
model data described above to produce the special mass property quantities 
required by the UFSSP for each flexible body of the spacecraft model. In 
particular, the following quantities are calculated: 


" J ■ / 

dm'' = Mass 

of Body j 


(9-1) 

•v 





dJ . i. 

f r'' dm'' 

■ t 4 4 - 

3=1 

- P 1 

(9-2) 

P = -L- 

* m'' 

f H dnrl 

B^ 

• Z *1 6 4 

3=1 

(£=1,2,3. . .n ) 

(9-3) 

^ = - 
1 p 

I H 

B J 

■ Z Y is 

3=1 

P (£.= 1,2 ,3. . .n ) 
-p J 

(9-4) 

z j = — 
k l j 

m 

/ +k * *1 

B d 

3 

dm'' = (k, £=1 , 2 , 3. . .n^ ) 

3=1 

(9-5) 

A 

*<=> 

ii 

3 -r 

/ 

dm^ (k, 1=1, 2 , 3. . 

■ n.) - t^/w? 

(9-6) 

b] = K 
1 P 

/ ♦] 
J j 
B J 

-z z 

a=l 3=1 

B Lb 4 sj 

(9-7) 


(£— 1 ,2,3,. . . n^ ) 


102 



(9-8) 


where 


where 



1 


nr 



dm j = 


t i 

u=l 0=1 


C j e j e j 
k£a8 -a -g 


(k, £=1,2,3, • • * n j ) 





3 3 

-£ £ 


a=l 8=1 


r^ • dm^ (£=1,2,3, .. .n.) 
5 (r^ • - r^ dm^ 

8 Ls 4 4 (H-1,2,3,...^) 


N' 


£ag 



B 


j 

£ag 




a=l a=l 


+i> - +| ‘lnr> 

e g (k , £=1 ,2,3, • • • n j ) 





C j 

k£ctg 


(9-9) 


(9-10) 


(9-11) 



(9-12) 
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where 1^ is the standard inertia tensor for Body j with respect to the 
e^ frame with origin at the Body j hinge. 

“P 


In several of the above equations , 


6 


3 3 


■ E E 5 «b 

o-l $-1 


e J 

-a 


4 


Kronecker delta tensor 




(9-13) 


Finally, the generalized stiffness matrix for Body j is computed 
as follows: 




j 5 k& 
m 


u 


U) 


(9-14) 


so that Kj^ is a diagonal matrix whose Jtth diagonal element is the square 
of the Ath eigenvalue for Body j divided by the mass of Body j . 


The integrations indicated in Equations (9-1) through (9-12) are 
replaced by a summation over the nodes (joints) of each flexible body when 
the indicated quantities are calculated. 
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X. 


Definition of Quantities 


The following tables define the quantities used in the UFSS Program. 
"Body (F,R)" specifies if the quantity is used for a flexible or a rigid 
body; "Comp" specifies the program area in which the quantity is computed 
and "Used" specifies those program areas where the quantity is used — 
according to the following code: 


1. 

Orbit Subroutine 

7. 

Flexible Combining Algorithm 

2. 

Disturbance Subroutine 

8. 

State Vector 

3. 

Control Subroutine 

9. 

Mass Properties Subroutine 

4. 

Auxiliary Calculations 

10. 

Input 

5. 

Sequencing Algorithm 

11. 

Initial Conditions 

6. 

Rigid Combining Algorithm 




The units are mass (M) , length (L) , time (T) and the derived unit force (F) 
-2 

equal to MLT 
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Appendix A » Derivation of the Basic Dynamic Equations 


This appendix contains the derivation of the basic dynamic equations 
for the UFSS Program. The first section gives details of the derivation 
of the governing equations for a single flexible body. Section A. 2 docu- 
ments the derivation of the Flexible Combining Algorithm, while Section A. 3 
documents the derivation of the Rigid Combining Algorithm. 

A. 1. Equations of Motion for a Single Flexible Body 

Since only terminal bodies of the spacecraft model are assumed to be 
flexible } the equations for an arbitrary flexible body (Body j) are 
derived explicitly considering its limb (Body i) . Specifically, referring 
to Figure A.l the instantaneous position vector a^ to an arbitrary mass 
element dm^ in Body j is as follows: 

= R r + R 1 + + P + u j (A-l) 

where R r is the inertial position vector to a reference axis frame, R 1 , 

, r^ and u^ are as given in Figure 2.3. 



Figure A.l. Position Vector Definition for Body j 
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The first time derivative of 


a ^ is 


dt k 1 


R r + R 1 + (m r x R 1 ) + Jl 1 ^ + (to 1 x l 1 ^) 


+ ( 10 ^ x r^ ) + + (to^ x ) 


(fr - 2 ) 


— r — — -1 

where to , to and to J are the angular velocity vectors of the reference 
r i 

axis frame (e^) , the Body i fixed frame (e^) and the Body j fixed frame 
(e^) , respectively. (The e^ frame is assumed to define the undeformed 
position of Body j.) 

The Kinetic energy of Body j, T^ , is formed as follows: 


1 f d d , -U 

' 2 J , dF ( “ > ' dF ( “ > 


dm J 


Therefore, substituting Equation (A-l) into (A-2) and performing the necessary 
integrations over Body j, the Kinetic energy becomes: 
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HZ 


scalar 


nr 
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) : 




second order unit tensor 


= Kronecker delta 


) yields a second order tensor. In particular, the 
following identity has been used in determining the above expression for 

T J . 

(A x B) • (C x 5) = (C • A) (5 • B) - (C • B) (5 * A) 

= C • [6(B ‘ 5) - 5 5} * A 

where 


A, B, C and D are vectors.] 

The strain energy of Body d is defined to be 

= |» J ii 4 4 4* 

k,£=l 


where 


= the generalized stiffness matrix of Body j associated with 
the generalized coordinates and q^ . 

Likewise, the dissipation function for Body j is defined to be 


n, 

3 


E 4 4 4 

k,Jt 3 l 


kJ i 
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where 


„ = the generalized damping matrix of Body j associated with the 

k& , 

generalized velocities q^ and q - ^ . 


of R 1 


Now taking as generalized coordinates for Body j the three components 

the p gimbal angles 0^ (y=l,.. .,p.) and the n deformation 

j J Y J J 

coordinages q? (k=l, . . . ,n . ) , the Lagrange Equations for Body j are as 
K j 

follows : 


d_ 

dt 


d_ 

dt 


d_ 

dt 


3T J 


3q: 




3T 


30 J 
• Y- 

3T J 


3R 1 
u a J 


5T- 


3T J 


30- 


3T 


3U J 

*4 


3U J 


30 J 


3U' 


3R 3R 

a a 


3D J 

3q i 


- <3.7 (k=l,...,n) 


3D 


30- 


3D J 

3R 1 


. - (y 53 !* • • • >P.s ) 


(a=l,2,3) 

a 


(A-3) 


where , © j and are the generalized forces (environmental distur- 

ic “Y ot 


bances, control forces, etc...) associated with the generalized coordinates 

qrj , and R^ , respectively, 

k y Y 

A- 1.1 Equation for q^ 


Performing the indicated partial differentiation, the equation for 


J 


q, becomes 
k 

j 


m 


+ 

t 1 

+ 

• , 

2(u r x R 1 ) 

+ (u> r 

x R*) + </ 

X 

(to 

S 1 * 

+ 

2 G 

^ x ) + 

(m 1 

x Ph 



-i 

0) 

X 

(m 1 

x J • 

4 + 

4 t 4 + 

• 

K 

S J - 

u 

q j 
- i 

- 2 • 

' z L 4 

- » J ■ 4 

• 

toj 


• 

Kt 

• “ J 4 + 

’ 

i + v m 4 1° 

£ 

j 

k 
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(Note in the above equation and in the following, repeated subscripts are 


summed over their allowable range; e.g.. 


4* 4 


2 ^ 4n4’> 

£-1 


Rearranging and writing the above in component form 


we have 


m 


*? A 

ka _ 


je 

ag 


4 + 


4s 4 + 


ag 


(2 u 


3p 


R 1 + 
P 


w 


0P 


- 4> 


r 

’ey 

~r 
oj . 
y6 

E s> 

+ 41 

~i 

“ey 

i 

ai jc 

y6 

+ 41 


+ 

2 *>* 
ey 

ka 

A jl 
' a6 

l ij 

sg 

“6 + 



+ yI b? - 

ka a 

ka 

a 





- 2 H^ 2 0)^ 

ka a 

- 

N kag ' 

+ q^ 

k£ q £ 

1 

= «3 ^ 

^k 

where : H.^ = 

ka 

q i 

Z j 

k£a 

II 

CVJ 3 

q i 

Z j 

k£a 

“kaB ■ 

4 

gj 

k£ag 


41*4 4 + 4 * 


q 


j 

i 


Rearranging once more. 



+ sj 4 ^ 

SL ka a 

- 

ka 


w* + 

6g 3 ka 

A j « i 

ag 3 



+ 

<t>^ 1 + 

ka a 

ka ag 

f 

- 2 Hjj 2 a) j - 

ka a 

S jl ° + 
b k + 

4t 4 

(A-4) 

+ 

v « 4 j 

= g> $ - 


4 £ 

ka ag g 





where: 
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s jl 


OC 

Ji 


< 2 ?g R g + %B R 0 + % V R 3 ) 


+ a j * oj 1 n* j 

aa a p pB 3 


i ^ 2 = + 2 & 

a a3 B 


o 

S j4 

kg 

.jlO 


= Y: 


kg 

j 


(N: 


"ii 


,35 


,j 


kag 


+ W < *4 


Now, introducing the following relationships between w 3 , oj 1 and the 
gimbal angle derivatives. 


where 


and 


J 

- A* 1 , 

HJ 1 + 

g{ 

8^ + 

*J Jo 

G^° 

§Jo 

8 

8a 

a 

By 

Y 

3<S 

5y 

Y 

'J 

= aJ 1 , 

i 1 + 

G o 

“i 
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6 

Ba 
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By 
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6 
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A 

JoiW 


8 
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ya 

By 

Y« 
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By 

Ya 
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Uo 


the equation for q, becomes: 

l & 
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K* 
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A j q « 
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J, 1 + 

$ 3 

A^ r R 1 

+ 
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8 

ka 

aB B 


ka 

+ 

ka 

s 3? 
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4 10 

+ 
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+ 

V 3 
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(A-5> 
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where 






Finally, the equation for 


can be written in matrix form as follows: 


3 


All A12 A13 r A14 

B kH B ky B kg B k$ 


] 









C 


A1 

k 


where 


r A11 

k£ 

= m^ 

"k* 


A12 

B kY 

= m J 

s ii 

ay 

r A13 

B kg 

= m j 

sf 

ka 

4 

ag ka ao 6g 

M4 

B kg 

- m j 

ka 
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ag 
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°k 

■ - “ 3 Hi 

S f + »ka S a ? 
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hP 2 
\s 

4 - 4 10 


+ 

4 * 4 

+ V k£ ’ll + «£“ 


It should be noted in the above equation that 



( . JP-.. 
\ CR r ) 3 


ka ag 


4 


(A-6) 
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that is, the gravitational force acting on Body j if Body j were located 
at the origin of the reference frame has been removed from the external 
generalized force term by equating it to the corresponding acceleration 
term to produce the following orbital equation 


Ym 

<R r > 3 


A J ! R r 
kct ct& 


. - «J «2 £ 

$ ka a3 $ 


or 


(A-7) 


B). 

forces associated 

Body j had its 

axis frame e r 
-a 

(terminus of R"). 

A. 1.2 Equation for 

Once again performing the indicated partial differentiation, the 


equation 

for 

R 1 

a 

is as 

follows: (note 

that 

R 1 

= 

R 1 e r ) 
a -a 

j r , 
m e 

-a 

’I* 

4- 

-i 

R 1 + 

2 (ui r 

x R*) 4 

• (w 1 

X 
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4- 
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X 

(m r x 
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-ii 

SL 2 4- 

2(m i 

X 
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4- (S 1 x I i:3 ) 4- 

+ 

-i 

0) 

X 

(oj 1 X 

• 

& ij ) 4- (m j x 
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+ 


x (bP x d^) 4- 

+ 

*k 

"j 

q k 

+ 2(w j 

x P k ) qj + 

<3 

X 

*k> 

4 * 


+ 




X 




* 


1 a 


* (R r ) 3 0 

which is solved in the Orbit Subroutine (See Appendix 


i 6 

Thus, represents all external generalized 
with q^ except the gravity force which would act if 
mass center positioned at the origin of the reference 


(A-8) 
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m J 


Rearranging and writing the above in component 

form we 

have 

R 1 
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* 2 “aB 

Ri 
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where H^ 3 = <l4 . Again introducing the relationships (A-5), the above 

p K Kp 

equation becomes 
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or more simply 
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where 


S j6 = + H^ 3 

a a a 


H 


j4 ■ *2 * 2 . 


Finally, the equation for can be written in matrix form as follows: 
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A31 „A32 „A33 „A34 


aSL 


where 
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ay “8 ag 


] 




Y 
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(A-9) 
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Once again, somewhat as in the case of the equation for q 
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Thus, the gravitational force acting on Body j if Body j were located at 
the origin of the reference frame has been removed from the external gener- 
alized force term by equating it to the corresponding acceleration term 
to produce the orbital equation 


or 


ym J 


(R r ) 3 


A re R r 
a$ K B 


*3 


= o 


~Y 


(R r ) 3 


= R, 


which is solved in the Orbit Subroutine. 


In addition, the hinge force acting on Body j due to Body i (F^ => 

F^ e^) has explicitly been separated from the remaining external gener- 
alized forces in order that it may be most conveniently utilized in the 
Flexible Combining Algorithm to be described shortly. 

A. 1.3 Equation for 

The equation for 8^ can likewise be obtained by performing the 
"appropriate" differentiation. However, in this case T^ is an implicit 
function of 0^ and 0]j only through the appearance of and various 

transformation matrices occurring in the expression for T^ when written 
In terms of components of the vector quantities. The required differentiation 
and subsequent algebra is extremely long and a much simpler formulation 
arises by considering the torque equation about the Body j hinge: 



+ 



x 





fj + ji 1 


(A-10) 


where T" 1 is the total external moment on Body j about the Body j hinge and 
i ji is the hinge torque acting on Body j due to Body i. The second derivative 
of is 
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d 2 -j 
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dt 2 


= 

r 

+ R 1 + 2 (m r 

X 

i 1 ) 

+ 

K 

u 

•1 3 

N— ' 

R 1 ) 


+ 

-r 

ai 

x (aj r x R^ - ) 
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•- • • 

+ 2(o)^ x u d ) + (oj^ x u^) + a ^ x (a>^ x u^) 


Substituting (A-ll) into (A-10) and performing the necessary integrations 
over Body j , the torque equation becomes 
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Introducing the relationships (A-5) , the equation for 6^ becomes 


m 


ai 


where 
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for 
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: written in 

matrix form 
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where 
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Here again. 
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j) - T ji . 

6 / a 


Thus, the gravitational force acting on Body j if Body j were located at the 
origin of the reference frame has been removed from the external generalized 
force term by equating it to the corresponding acceleration term to produce 
the orbital equation 


Ym 


j 


(R r ) 3 


3 J6 A j ® R r 
aB 86 6 


R r 

aB A B6 k 6 


s j r AiT r; - m j s j : mz r; = o 


or 


Rt = - ~ m -~ R 




CR r ) 3 


which is solved in the Orbit Subroutine. 

_ j i 

In addition, the hinge torque acting on Body j do to Body i (T J » 

T^ 1 e^ ) has explicitly been separated from the remaining torques in order 
a -a 

that it may be most conveniently utilized in the Flexible Combining Algorithm 
to be described shortly. 
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A .1.4 Constrained Degrees of Freedom 


1^ ID AM J 


The equations just presented for , R^ and 0^ have been derived 
assuming that all three degrees of rotational freedom exist at the gimbal 
hinge joining Body j and its limb, Body i. In general, there will exist 
fewer than three unrestrained gimbal axes and the equations must reflect 
this fact. 


Since the transformation matrix is determined from three sequential 

Euler rotations through the angles 0^ , 0^ and finally 0^ (see Section 4.1), 

the order in which the gimbal angle rotations are locked is important. 

Thus, 6^ is constrained first, 0^ is constrained nect and finally, if 
necessary, 0^ is also constrained . Thus, if a single degree of freedom 

exists between Body i and Body j, it is 0^ about the first gimbal axis 

j q i 1 

e^ ; if two degrees of freedom exist, they must be 0^ and 0^ • 

In the equations (A-6) and (A-9) for q^ and R 1 respectively, 

constrained degrees of freedom are simply handled by replacing by 

j-fr* 

G wherever it appears in the equations, where the super + implies 
ay . 

removal of the Xth column if 0^ = 0 . 

In case of the equations for the 0^ themselves, constraints are handled 
i "h 1 Y 

by introduction of G wherever G apears in addition to elimination 

ay ay 

of the equations for the constrained variables. This latter operation is 
accomplished symbolically by introduction of a super-zero notation as 
follows : 
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F-re: <$» Y = l,--., Pj ; o, 0 = 1, 2, 3; k, 4 - 1 n 

where is the number of flexible degrees of freedom of Body j and 

is he number of relative rotational degrees of freedom of Body j. In 


J 



case p. * 0 , the above equation is void and io^ is determined directly 

-i J J 

by to and the input values of the 6^° through equation (4-10): 



Details of the Flexible Combining Algorithm 


The Flexible Combining Algorithm is used only when System A in a 
given combining operation is a single flexible body, call it Body j. 
Referring to Figure 5.1 > assume that System B has Body i as its member of 
lowest level (i<j) and that System A is to be connected to System B to 
yield System C. (See Sections II and III for a detailed description of the 
system model and notation.) Assume also that Body l is the lowest numbered 
branch of Body i in System B and that j<£. 


The equations for System A are given by (A~6) , (A -9) and (a- 12) which 
are written below for reference 



Let us now assume that the System B equations are of the form 


(A -13) 
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System B 
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b B33 

aB 
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c B3 

a 


w i£ 

a 


i £ 

where L is the total hinge torques on Body^i exclusive of the torques 

01 i £ 

imposed by its branches numbered i £ and W is the total hinge forces 
on Body i exclusive of the forces imposed by its branches numbered k £ . 

Here, B „ , B . , B „ and B . are 3x3 matrices with the 
oiB aB aB aB B2 B3 

indices a and B running from 1 to 3; C and C a are 3x1 matrices; 

"it , „ -Bll . “ d B12 , B13 

if 0 has r components, then B is r x r , B „ and B „ are 
n g21 B31 mn B1 mB mB 

r x 3 , B and B are 3 x r while C is r x 1 . As used above, 
an an m * 

0 n contains a component for every flexible and rotational degree of freedom 
of Body i's branches and sub-branches numbered greater than or equal to £ . 


(A -14) 


We shall now show that the combined System C equations have a form 
identical to the System B equations, namely 


System C 
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c c3 
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a 


(A-15) 


and exactly how the elements of the coefficient matrices of (A-15) are 
synthesized from the coefficient matrices of (A-13) and (A-14) . Having 
accomplished this, it will be obvious that the same synthesizing procedure can be 
used to combine an additional flexible body to Body i since the new System B composed 
of the original System B of Figure s.^plus Body j is governed by the Eqs.(A-15) 
which are identical in form to those of the original System B as given by (A-14). 
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The required equations governing the motion of the combined System C 


are as follows: 

A) equations of motion for Body i translations (A-16) 

B) equations of motion for Body i rotations (A-17) 

C) equations of motion for Body j along those gimbal axes where 

it possesses a degree of freedom with respect to Body i (A-18) 

D) equations of motion for the branches and sub-branches of Body i (A-19) 

E) equations of motion for the flexible coordinates of Body j (A-20) 

Let us first consider the equations for Body i translations. The 
hinge force relation for Body i is given by 


w iS ' = w 1 ^ + F^ . 

But, from Newton's Third Law 



Therefore, using (A-13) to eliminate the forces of interaction we have 
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Thus, utilizing the last equation set of (A-14) 
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or finally 
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Thus, the following elements of (A-15) are now determined: 
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C32 

B32 

B ag 

+ A ir 

A33 

B ag 

ag 

i A 

a a 

C33 

r B33 

B ag 

+ A ir 

_.A34 

B ag 

ag 

a a 

,C3 

a 

c B3 

a 

+ A± « 
ag 

r A3 

C g 



where, necessarily: 



"i£, 

0 

n 




e J 


with 


(A-21) 


(A-22) 


a, 6 = 1,2,3 

b = 1,2 (r + n j + Pj) 


and the column matrix O' 1 

Y 

the hinge interconnection 


containing only the degree-of-f reedom components at 
of Bodies i and j. 
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The equations of motion for the branches of Body i are exactly those 

“i £ 

given by the equations for 0 in (A-14) : 

b B11 + b b j 2 ;; + b B33 % - c B1 . 

mn n mg $ mg g m 

The equations of motion for the flexible coordinates of Body j are 
exactly those given by the equations for q-J in (A-13) : 

All "j A12 "j A13 ;i A14 R i _ Al 

B k£ q * + B ky \ + B kg “g + B kg R g " C k • 

The equations of motion for Body j along those axes where it possesses 

a degree of freedom with respect to Body i are exactly those given by the 

equations for 6^ in (a~13) : 


A21° "j A22 

B 6Jl q * + B 6 Y 


P + b a23 


•i + r A 24° "i _ „A2‘ 

“g + B 6g R g ‘ C g 


+ T 


ji 


Y <50 0 6$ “g 6 

Therefore, the following elements of (A-15) are now also determined: 

Bll 


B 


Cll 

ab 


B 

mn 


B 


0 

0 

,A11 

b£ 12 

ky 

ki 

A21° 

-A22 

B «Y 

'61 


C12 

ag 


B 


b: 


B 


B12 

mg 

A13 

kg 

A23 c 

<$$ 





B!3 

mg 

1 

„C13 

B ag 

h 

r M4 

B kg 




r A24° 

B <Sg 




J 


(A-23) 


(A-24) 


(A-25) 
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<? 

of * 


Finally, consider the equations for Body i rotations. The hinge 
torque relation for Body i is given by 


L iA = L 1 ^ + T 1 ^ + x 

But, as before 



and similarly 

_ _ ^j i 


Thus, from (A-13): 


F iJ - - a 1 !; 

a ag 


_ A3 A31 ..j 

C g + B 0i q * 


Bi 32 » + B A33 i* + B A34 


ftY 


'ftp 


’ftp 


and from (a- 12): [note that all three rotational equations for Body j are 

required here to eliminate the interacting torques] 


ij _ ij _ C A2 A21 "j A22 gj A23 ‘i A24 

r a a ub c g + b 6a q a + B av y ap % + B a P 


Therefore , 




ij 


A ij 0 

! - C A2 

+ b a21 

ag j 

i 6 

a* 


^ !- 

cf + 

aa 

<*3 

a 


By 


A23 -i A24 gi 

B ep % + B ap R p 


,A31 "j + b A 32 gj + b A 33 -i + B A34 ”i 

ea $y y ep p ap p 
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Thus, using (A-14), the equations governing Body i rotations are: 


b B21 e u 

an n 

+ 

n B22 .i 
B «6 “8 

+ B 823 L 1 = C B2 + L 1 * 
ag g a a 

+ i 

ag / 

A2 _ An “j _ 
C g B g£ q A 

Bf 2 V - B A23 m 1 - 
By Y Bp p 

+ & A^ic 

aa ag j 

f - « 

A31 

g£ q £ 

- C i - i 33 i ■ 

or, more simply: 



b b21 5 1 * 

an n 

+ 

it ’i 

+ o 2 

•ay 

i + it " it i 

= Q 5 + 
a 

L 1 * 

a 




where : 





it 

= 

a ji ( t) 

ag 

B A21 + 
g£ 

oij A lr E A3I 
a a ag U 

Q 2 

ay 

= 

A j i (T) 

ag 

sf 2 + 

By 

£ 1J A 1 * B A32 

aa a g 3y 

Q 3 

^ag 

= 

b b22 + 

ag 

a^< t > 

aa 

B A23 + A lr B A33 

ag aa ap pg 


= 

B B f + 

ag 

A ii(T) 
a a 

B A24 + A ir B A34 

ag aa ap pg 

< 

= 

c B2 + 

a 

A^ (I) 

a8 

C A2 + A 1 * cj 3 

g aa ag g 

Thus , the 

remaining elements of (A-15) are now determined 

C 21 

B ab 

= 

[ b b21 

|_ an 

t I 

I o 1 ] 
1 '<** \ 

q »y] 

C22 

B ag 

= 

i e ■’ 

R C23 

B ag 

’ Q »S 

c c2 

a 

= 

i 




B 


Bp P 
A34 -*iJ 


- B 


Bp 


(A-26) 


137 



Having now demonstrated that, given the System B equations in the 
form (A - 14) , the combined System C equations are indeed in the similar 
form (A-15) , it remains to prove that the System B equations do indeed have 
the form of (A-14). There are two cases to consider: first. System B is 

a single rigid body and; second, System B is an arbitrarily interconnected 
system of rigid and flexible bodies with Body i being rigid. 


Considering the first case, one simply has an initialization of the 
System B equations. Here, the appropriate equations are the familiar 
Newton-Euler Equations : 


m 1 jl r 

+ 

R 1 + 2 

• . • 

—r — x 

co x R + ui x 

=i 

-i 

-i 

=i -i -i 

I 

0) 

+ 0) X 

I r ui = T 

where 

i 

m 

and 

are respectively the 


-i 


-r 


,-r -i, 


F 1 


(A-27) 


(dyadic) for Body i while F 1 and T* are respectively the total external 
force and centroidal moment acting on Body i. 


Resolving the above equations into component form and eliminating 
the orbital equation from the force equation, the Body i (System B equations) 
become : 


m 


A ir 

u0 

l®e + 

o r 

2 a) „ 

00 



r 

= 

F 16 + 

a 

w ij ■ 

a 

•i 

“0 

+ 

ab 

T 1 

1- - ai 

B6 


*6 + 


-v F 


ij 




06 


u> 


T le + L ij 
o a 


By 


r . 
a) ,.) 
yS 


riJ 


ag *8 


(A-28) 


where F is the total external force acting on Body i exclusive of all 
a 

hinge forces and the gravitational force Body i would experience if it 

-j- 

were located at the terminus of R , while is the total external 

centroidal moment acting on Body i exclusive of all hinge-produced moments. 
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Thus, in this first case the System B specification [i.e., specification 
of the coefficient matrices of (A“14)] is as follows: 


B 


B 


B 


B22 

aB 

B33 

«3 

B23 

<*B 


a0 


j A ir 
nr A . 
aB 


= B 


All the remaining sub-matrices of B 
components (m,n=o). 

In addition, 


A3 3 

aB 

B 


- 0 


"i£. 

are void since 0 


has zero 


„B2 


,B3 


- u 


«B 


= - m 


i« 4 + T le 

By <$ a 


^ ir 

aB 


U) 


36 


4 J + 


while C B ^ is void, 
m 


, *r 
(a) B6 


r r s 

111. w - J 

BY y6 


+ F 


le 


Considering now the second possible case; that is that System B is 
an arbitrarily interconnected system of rigid and flexible bodies with 
Body i being rigid, it is sufficient to demonstrate a second combining 
algorithm which synthesizes the equations for two arbitrary systems of 
interconnected bodies when System A is not a single flexible body. This 
algorithm will henceforth be denoted as the "Rigid Combining Algorithm", 
and its output will now be shown to have the required form of (A-15). 

A. 3 Details of the Rigid Combining Algorithm 

Consider the combining of two systems (Systems A and B) to form a 
third system (System C) as shown in Figure 3.3. Here, Body j of level 
(N+l) is the lowest leveled body of System A and Body s is the lowest 
numbered branch of Body j. Body i of level N is the lowest leveled body 
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of System B and Body A is the lowest numbered branch of Body i. Since 
the highest numbered branches of a given limb are connected first by the 
Sequencing Algorithm specification, it is necessarily true that j<A 

Here, System B is with no loss of generality identical to the 
System B of Figure (A-2) , hence its equations must be identical to those 
given by (A*14), namely 


System B 






“ A 


— — 


— - 

b b11 

mn 

_ bi2 

mS 

R B13 

B mB 


e i£ 

n 


c B1 

m 


0 

b b21 

an 

— B22 

b «b 

r B23 

B aB 


•i 

tt B 

n 

c B2 

a 

+ 

A 

L 1 * 

a 

b b31 

an 

r B32 

B aB 

r B33 

B aB 


R 1 

r b 


c B3 

a 


a 



— 


— - 


— — 


— — 


In addition. System A must have a similar description and interpretation 


(A- 29) 


System A 


All _A12 A13 

\a \b % 


A21 A22 A23 

V B aS uS 


A31 R A32 A33 
B aA B a S B a $ 



V 




0 



= 

c A2 

a 

+ 

A 

L Js 

a 




c A3 

a 


A 

a 


— — 


— — 


- - 


Here, B 


A22 




-A3 2 , _A33 

B „ and B 


'aB anc * ~a3 
M components, .then 
,A31 


and B 


aA 


are 3x3 matrices; 
All 
kA 

are 3 x M while 


Bf 


is M x M, 
„A1 


,A2 

a 

a 


J k 


is M x 1. 


A23 

A3 "is 0 ® 

C A are 3x1; if has 

a A13 A21 

and b£“ are ,M x 3, B^ 1 

Once again, 0^ 3 contains a component for every flexible and rotational 
degree of freedom of Body j’s branches and sub-branches numbered greater 
than or equal to s . 


and 

A12 

kB 


(A-30) 
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We will now show that the combined System C equations have the form 
of (A-15) as desired. Specifically, the required equations governing the 
motion of the combined System C are as follows: 

A) equations of motion for Body i translations 

B) equations of motion for Body i rotations 

C) equations of motion for the branches and sub-branches of Body i 

D) equations of motion for Body j along those gimbal axes where it 

possesses a degree of freedom with respect to Body i 

E) equations of motion for the branches and sub-branches of Body j 

Before deriving the above equations, certain relationships must be 
established. Since Body j is a rigid body it is true by definition that 


R j « R 1 + I i;1 - P 1 . 
Thus, taking a first time derivative, 

P + ui r x R^ = R 1 + <j r x R 1 + + w 1 x 


- P x » Ji 


or 


R J = R 1 + a) r x (R 1 - R^) + + co 1 x i, 1 ^ - ^ x . 

In component form (and substituting for R^ - R^ ) one obtains the recursive 

• i • i 

relationship for R*' in terms of R : 


R J 


+ A 


- A 


+ A 

ir(T) 

aB 

•ij _ jr(T) 
*B aB 

P 1 

B 

ir(T) 

ctB 

"i 

U) 

»Y 

~r .ir(T) ( 

~ “aB A BY j 

: & 
Y 

jr (T) 
oB 

“by 

~r jr(T) | 

aB A By ( 

Y 


(A-31) 
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Taking a second time derivative in vector form: 


1* 

+ 

-r 

to 

X 

-i 

R J = R + 

-r 

to X 

R 1 + co 

r x (E 1 

- R j ) + u r 

x (R 1 

- P) 


+ 

-r 

to 

X 


x (R 1 - 

- ) 

] + 

^ -i 
+ to X 

pi + to 1 x 

pl 



+ 

-i 

to 

X 

~pi 

. -i 

+ to X 

(2 

x I ij ) 

- P 1 - 

-J x i ji - 

P X 

P 1 


- 

5* 

X 

P 1 

- X 


X £ jl ) 





or. 

in 

component form: 









"1 

R J 

a 

= 

"i 

R + 
a 

_10 -i 
P «6 “6 

+ P 

19 ^ 

. (0„ + 

ctS B 

P 16 

a 



(A-32) 


where 


,10 

ctS 

,19 

ot$ 

-.16 


-A ir(T) 

ay Y0 

A Jr(T) 


ay 

aS 

aS 


V 8 


P* 

S 

+ 2 »U 

PI 

6 

+ 

2 

By 

?! 

V 1 

B 

+ 2 “e« 

P 1 

+ 

0i4 

By 

?| 


•r ~r ~r / Li 

to . + to to „ > <R„ — 

aS aY YB [ j B 


'8 


+ 2 


'a$ R S 


- R' 


Let us first consider the equations for Body i translations, 
hinge force relation for Body i is once again given by 


The 



ijlJ + plJ 


W 


so that, using the relation 




, one has 




A j ; (T) 4 1 * 

a$ B 
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But, since Body s is the lowest numbered branch of Body i,it follows that 


W jS = F^ 
a a 


so that 


W 


iJl 


W 


ij 


- A ji(T) W jS 
ag W B 


(A-33) 


Thus, using (A-29) and (A- 30), 


B31 

'an- 

e 1A + 

n 

<? 

•i ^ 

U)g + 

R B33 k 1 
B ag R g 

C B3 

a 

ji(T) 

ag 

Hi 1 


+ b a32 
+ B gs 

*j + b A33 
“6 + B g<S 

pj 

6 


i . 

Substituting for and R“ one finds (noting that here P 

^ + >1? “l ♦ »S* 3 - 'I' - ? 

- A »e <T> | “el 1 + + «£ *1 + p s 7 ] 

[*i + + (< + < *1 




w 


ij 


+ B 


A3 3 
B6 


S^ 5 ) 


* ■;')* >.“KI 


or 


an n 

B32 
ag 

B33 


+ [b 

+ [ B 


+ A' 


ji ( T) 

ag 

.A31 

V 

B i 9 

+ a^ ( t) 
ag 

fB A32 

(B g6 

G i + 
G 6Y 

+ 


p ” 

G* + ) V 
PY Y 

ji(T) / 
“a6 \ 

.A3 2 


+ b A33 p“ 
6y YB 

+ 

r A33 

B 5Y 

P 19 

YP 

A* 1 ' 

P3 y 

)]"6 



ji(T) A33 
a 6 B 6g 


]*B ’ + 

« 

A33 / 19 
+ B 66 V 6P P 


C - 

a 

A J 

ag 

17 + 

p“' 


B A32 P 17 
gS o 


or 
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b b31 

an 

"1Z 

0 + i 

n 


A31 

8i 

®i s + 

p 12 

a8 

+ P 12 

a<$ 

Ji ( *i 
A S6 (“6 

♦ ! «u 3 

+ 

A ji(T) 

06 

A33 
B <5 0 

ji(T) j A3 
ag | C 8 

t>A33 
" B 86 

p 16 1 

6 i 

1 ±j 

► + W J 
1 a 


where 






P a6 = 

A jl(T) 1 

a 8 j 

1 A32 
l B 86 + 

A33 

B 0Y 

p “! 



J+ 

'By 

eJ 

Y 

+ 1 

B + 

otB 

A ji(T) 
a S 

R 1 


c B3 

- P 12 

P 17 

R 8 


a 

P a0 

P 8 


B 6y P Y 8 


(A 


Finally then, the following elements of (A-15) are determined: 


C31 

ab 

- b b31 

1 an 

! 

i 

i 

A jl(T) 
a 8 

A31 

B 82. 

C32 

a0 

- r B32 

' B a0 

+ 

JUT) 

aS 

r A33 

®«Y 

C33 

a0 

- r b33 
■ B a0 

+ 

A ii (T) 

aS 

r A33 

®68 

C3 

a 

„B3 
= C 

a 

- 

P 12 P 17 

a8 8 

+ . 


I P 1 ? G^ + 


I 


ctB BY 


p 1 ? * pl * 

Y$ a6 60 


Ji(T) 

a$ 


A3 W A33 16 

C 0 ' V p 6 


The equations for the branches ^and sub-branches of Body i are exactly 

" 1 2 , 

those given by the equations for 8 J in (A-29) : 


b b11 e il + b b J 2 + b b J 3 rJ = c B1 

mn n m0 6 m0 0 m 


(A 


The equations for the branches and sub-branches of Body j are determined 

“■j g 

from those for 0 J in (A-30) : 


All “Js r A 12 -j A13 "j 

B k£ 9 Jt + B k0 “0 + B k0 R 8 


= c: 


A1 


* i ** 1 

Substituting for o)q and one finds 

p p 


.-34) 


-35) 
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All ys ’ ,A12 ji ;l j+ "j .17 

B kl 6 1 + A e« “S + G 6Y \ 6 


- C K - ^ 


( . ji *i 
A J co 
ycr a 


+ g^ + 0^ + P 17 

ya a y 


) 


+ P 


16 


= c: 


A1 


or 


B *11 ;j. + pi jj + P * ;1 + B*i 3 l 1 . 

ki i ky y kg P kg B 


= P 


13 


(A-36) 


where 


ky 

kg 

r.13 


R A12 + R 

\g + B 
A12 


c f ' j \s 


+ B 
A12 


A13 

P 19 ! G^ + 

k6 

6g f b gy 

,A13 

P 19 !a J1 

ko 

P crS | ag 


+ B 


A13 „10 


k6 


5g 

A13 „16 


+ B ks 3 p fg p g 7 - B kr p 6 


The equations of motion for Body j along chose gimbal axes where it 
possesses a degree of freedom with respect to Body i are obtained from the 
equations for in (A-30): 

p 


„AU ys + „A22 -j + b A 23 Jj 

aSL Si ag g ag g 


C B2 + L^ 3 

a a 


But, 


A 

L jS 

a 


T ^ 1 + F ^ 1 

a ag B 


and as above 


so that 


^ = W jS 

■ rs '' rt 


L jS = 1 + w| S 

a a ag g 


(A-37) 
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Substituting from (a - 30) : 


or 


R A21 

B ai 

e^ s + 
i 

A22 

B o8 

T ji 

a 

+ £ ji | 
ao / 

r A31 

B OJl 

A21 

B a£ 

- ^ 
ao 

R A31 

B oi 

r A23 

B a8 

- 

ao 

R A33 

B o8 


A + B A2 fl 3 2 
8 <*8 8 


- C 


A2 


+ b a32 a + b a p a 


08 


o8 


- C 


A22 

a6 


b A32 A 
ao 06 8 


j B A23 _ ~ji B A33j"j = T ji _ ‘ji C A3 

/ a8 ao o8 ) 8 ct ao o 


+ C 


B3 


A2 


• -j ** -1 

Substituting for and , one obtains 

p P 


B An . -Ji B A 3 lj 9 'jS + 
al a a csb \ l / a8 


A22 _ Jji B A32 

ao 08 


+ G 


+ P 


r e'J 

8y y 

19 
6y 


( a3± 

y YO 



- 8 31 C A3 + C 

ao o a 


A2 


- 2 1 b A33 

ap po 


11 * 


.ji *i 
a) 

I 8Y Y 


, p 10 -i 

+ P , a) 
6 y Y 


) 


a y ) + p i6 

6 


_ T 


ji 


Rearranging, 


p 2 e Js + p 3 
afl. I aY 


e 3 + P 5 + P 6 R a 

Y a8 8 a8 8 


= P 


14 


(A-38) 


where : 


2 

- Jtf 1 

- 2 

at 

c 

3 

. p 2 » 

G a + | 

aY 

| a8 

8y j 


P 20 - 

TJ A22 

B ae 


P a8 


A31 

o£ 


+ B k f pj» - ^ 

a6 68 ao 


A32 A33 19 

B 08 + V p 8 
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Where the super zero indicates retention of only the rows for which 0^ % 0 ; 
i.e., the Xth row is eliminated if 6^ = 0 . Thus, using (A-35) , (A-36) 
and (a- 38) , the following elements of (A-15) are also determined: 
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Finally, consider the equations for Body i rotations. Once again, 

+ + x F^ 


L U 


But, 


Therefore, 


T 1 ^ = - 1 and F 1 ^ = - F^ 1 


L 1£ = L 1 ^ - T jl - I 1 - 5 x F Ji 


and using the relationships 


1 = W jS 


T^ 1 . £ js - P 1 x F ji 


one has 


L 12, = L 1 ^ 


L* 8 + P 1 x W js 


£ iJ x W js 


Substituting from (A-29) and (A- 30) , 


b B 21 Jil + b B 22 ;l + b B 23 JJi . .B2 

an n aB 8 aB 6 a 


ij 


- A 


aa 


aa 


a l 


aa ap 


, A 
< s 

+ B A f 
cr B 

“b 

A23 

+ B „ 
aB 

"1 

r b 

- Cl 


A31 

pJL 

A 

+ 


% 

+ 

_A33 

B P6 

R i - 

c A3 

p 

> A f 

pi 


+ 


“ 3 B 

+ 

„A33 

B pB 

H - 

?\ 


or, 


(A- 39) 
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b e21 e 11 + P 18 B A f + A 31<1) B*? 1 

an n ap pZ ap pZ \ Z 


+ ^ap 8 B 


A32 + A ji(T) b a22 + 
pg ap pg g 


B B22 A + P 18 B A f 

ag g / ap pg 


^ji(T) b A 23 
ap pg 


+ B B23 2 - L 1J + C B2 + p 18 C A3 + A 33 '” cf 
ag g a a ap p acf P 


where 


“ 3 S 


B» 21 Jit + p? n. 
an n aZ Z 


18 

ag 

£ i;3 

ay 

A ji(T) - 
Y8 

A^ i(T > 

ay 

i ji 

Y3 

“1 

and Rp one 
3 

has 



jp 18 

/ ap 

„A32 

B P« 

+ a 31 < t > 

ap 

B A22 | ! 

P3 [ | 

i A ji *i 

A„ to 

I 3y Y 


+ G^ + 0^ + P 17 

Y 3 


+ P 


= L 


3Y 

10 

3y 

A 

ij 


+ B 


B22 

ag 


•i . p 19 

to + P 

Y 3y 


+ C 


B2 


+ P 


“b + 


p 18 b A 33 + A jl<X) b A 23 
ap pg ap pg 


R” 


( Ji *i 
A J to 
Y a a 


+ G J + y + p 17 ) + p^ 6 [ + b B23 £ 
Yd o Y / 3 ) Y3 3 

+ a 31(t) C A2 
ad a 


18 „A3 

«P °P 


where 


P' = P 18 B A31 + B A21 


,7 

a Z 


ap pZ 


ap 


P Z 


Simplifying the above. 


b b21 e 14, + P 7 8 jS + P 

an n a Z Z ay y 


8 + Ip 11 A ji n - P 9 P 10 


' ay yS 


aY *Y3 


n B22 -i , r B23 

% ( 6 + PoB 


9 "i 
P ag R g 


C B2 + P 9 P* 6 
a ag g 


p*: p 37 + a^< t > c A2 + p 18 c A3 + l 1 ^ 


,1! 

ag ‘g 


ag 


ag g 


where 


‘g 


(A-40) 


149 



A ji(T) b A 23 
aB BY 


P 

ay 






B 


A3 3 

BY 



P 18 B A32 + a^ 1(t) B A22 - P 9 P 19 
ay yB ay yB ay yB 


Thus, the final elements of (A - 15) are now determined: 


C21 

aq 

= [b b21 ! 

L an i 

1 > 8 1 
ai, j ay J 




C22 

'aB 

= - p 9 

ay 

P 1 * + P 11 

YB ay 

yB 

+ 

B22 

B aB 

C23 

aB 

- R 323 

_ B aB 

- P 9 
aB 




C2 

a 

= c B2 + 

a 

P 9 ^ ~ 

aB B 

P 11 

ctB 


+ A 3i ft (T) 
aB 


>18 

aB 



It has now been demonstrated that given the two arbitrary systems 
A and B of Figure C.3 with equations of the form (A-29) and (A-30), the 
equations of the combined System C are indeed of the form (A~15) which allows 
the combined system to be used as either a System A or a System B in a 
subsequent combining operation. Finally, if System A or System B consist 
of a single rigid body, the initialization procedure is identical to that 
given by (A~28) and the relations immediately following these equations. 
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Appendix B . Derivation of Environmental Disturbances 


B. 1 Derivation of Gravity Nominal Forces and Torques 

For a spacecraft in a Kepler orbit, the magnitude of the gravitational 
force at any point varies inversely as the square of the distance to the orbited 
body's center. Referring to Figure 2.3, one finds that the instantaneous 


force acting on a 
is given by 

differential 

mass element 

of 

the 

terminal flexible body 



d P 

Y P J 
(P j ) 3 

dm^ 



(B-l) 

where y is the 

orbited body' 

s gravitational 

constant , 



P j = 

R r + R^ 

9 



(B-2) 


~4 = 

R 1 + I ij 

+ 


+ u j * 

(B— 3) 

and 

(p j ) 

= (p j ’ P" 5 ) 

1/2 





13 -1 -i 

In order to determine (p J ) , one first calculates p J • p J as 


p j • p j = (R r ) 2 + 2 R r • R^ + (R^) 2 


so that 


P J * P J 


(R r ) 2 


1+2 


R r • R j f 

(R r ) 2 


+ 


R' 


r 2 

r 


and 


(pV 


(R r ) 3 


R r • R J ' f 2 


1 + 2 


(R r ) 2 


- 3/2 


(B-4) 
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A nominal value of (R^ / R r ) is 10 ^ so that (R^ / R ) is small 
compared to the first two terms in (B-4) and can be neglected. Expanding 
the resulting expression in a binominal series and retaining only first 
order terms yields 


<pV 


<R r ) 3 


1-3 


R r • rJ 

(R r ) 2 


Substituting (B-5) into (B-l), one obtains the following expression for the 
differential force acting on dm^ : 

— r — 1 

d F j - - — *-r 1-3 — + dl ” j ‘ 

(r v l (R > J 

Substituting the proper expression for r 4 from(B-l) one has 


d F J = 


(R r ) 3 


(R r ) 2 


Expanding the above expression. 


d F^ = - -w dm^ - ^ » (R 1 + + r^ + u^ ) dm^ 

(R r ) 3 (R r ) 3 


fiL. s r 

v L 


S r • (S 1 + + P 


(R r + R 1 + + p + ub dnA . 


Neglecting flexible displacements compared to the dimensions of the spacecraft 

IT -3 

and neglecting terms of order higher than (R ) , one obtains the final 

differential gravity force to be used in subsequent calculations: 


(R r ) 3 


— r- (R 1 + + r^ + u^) dm^ 

(R r ) 


dm^ 
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If Body j is rigid, then = 0 and the total gravity force is 


P = - V + F jG 

(R 1 ) 3 

where the first term is the force Body j would experience if it were located 
at the terminous of R r ; this term is combined with the proper acceleration 
term from the dynamics equation for Body j and solved in the Orbital Sub- 
routine as described previously in Section IV. Thus, 


F jG 


•/ 


d F j 


= - y m 


R J 


(R r ) 3 


+ 3 - ’ L ~ - (R r • R J ) R r 

(R r ) J 


But , 


R r = R r a r (see Orbit Subroutine) 


Therefore 


fJ g 


= - Y m ' 


j j 3 


(R 1 ) 3 


+ iJLEl . rJ) a r 
(R ) 


or 


r jG . _ 1 


Y 

r A jr 

R-^ - 3 

A j ! a* ( R j(T) A^ a^ )l 

(B-8) 

<R r > 3 

L ctB 

8 

cxB 8 \ Y Y8 8 )\ 


torque 

on the 

rigid Body 

j can be expressed as 



*jc _ 

f C 3 x 

J i 

B J 

d P 

(B-9) 


where p is the position vector of an arbitrary mass point in Body j. 
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Substituting the gravitational force in the form 


d F J = 


< 5 r + R - 1 + c J > 


(E r ) 3 


dnr 


+ 


3 Y 

(R r ) 5 


y • 


(R- 3 + c 3 ) 1 R r dm^ 




into Equation (119), 




or 


f (S' 3 x R r ) 
) _1 


(R r • dm j 


(R“)' 

B J 


f) G = i-JL- 
(R r ) 5 


f (R r • C 3 ) dnr 3 x R r 

J J 


B J 


Subtracting zero from the above equation. 


X3 G . -1 


f 

< E ‘> B J 

_ f 

(R ^> 5 y B j 

or, with I = the unit dyadic. 


It 


. R r j dm i x R r 


C 32 R r dnr 3 x R r 


T JG - R r x f {I (C j • C j ) 

(R ) 


R r - G 3 * R r ) dm j 


H-r R r x 

(R r ) 5 


/ 


B J 


{I (c 3 • i^ 3 ) " V S' 3 dm- 3 } * R 


rr 


3 - Y - R r X I 3 ' R r 

(R r ) 5 


where 1^ is the centroidal inertia dyadic of Body j . 
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Finally , 



_3JL_ 

(R r ) 3 


A 3 ! a r 


ctg 


A je 
05 5y 


4„<T) 4 


J e 


y a a p 


r 

a 

P 


where 
Body j 


I*' is the Body j inertia matrix with components expressed in 
ag ' “ 

axes with origin at the center of mass. 


B. 2 Derivation of Gravity Generalized Forces 

If Body j is a flexible body, then the general expression for the 
force on an element of mass, dm J , is 


d F 3 = - 1 — !_ dm d _ — 1 — (r 1 + I* 3 + r 3 + u 3 ) dm 3 
(R r ) 3 (R r ) 3 

+ 3 Y - - - {a r * (R 1 + I* 3 + r 3 + u^)} a r dm^ 

(R r ) 3 “ 

Therefore, the virtual work done by dF 3 acting on dm^ due to the virtual 
displacement 


is 


where 


also , 


6R 1 


6u“ 


6e J 


(? 3 


u j ) 




5u 3 


59 3 


Cr^ 




5R 1 = 

5R 1 

r 

e 


a 

-a 

5u 3 = 


*L 

59^ = 

69^ 

B ie 

a 

-a 


(B-10) 


(B-ll) 


(B-12) 
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<sw 


jg = R jg ^ + qjg 5| J + q(n)js 

a a x k ^k sv ~ / a a 


so that 


>js = 


/ 


d F 3 


r 

e 

-a 


Jg 


■/, 


d F J 


e 3 

-a ka 


Q(^ 8 = f** * { -i S x ^ + 

where R Jg ’ Q^ g and Q®^ 8 are the gravity generalized forces associated 

a fc a i 1 i 

with the generalized coordinates R^ , q R and £r respectively. 

Thus , 


R Jg = _ V m 

a 


(R r ) 3 « 


e r - ^-&r (5 1 + X iJ 

~° t (R r ) 3 


+ d 3 + H 33 ) 


r 

e 

-a 


3 vm J , r 
r3 (5 

(R ) 


(R 1 + + d 3 + H 3 3 ) }a r • e r 


-a 


or 


R 


jg 


_ JL 


m .re „r 


(R r ) 3 


» _ M . 

ag 8 


nv 


j 


(R ) 


{R 3 + A^ (T) H^ 3 } 

i\3 a ag g 


3y .re r 
(R r ) 3 


{(R^ (T) A re + h^ 3< ' T) A je ) a r } 


op 


Op' p 


Finally, subtracting off the orbital term as described in Appendix C, 

R jG = _ X_sl {R j + a^ (T) H^ 3 } 

a ^ R r^3 a ag g 

+ A re r {(R j(T) A re + R j3 (T) A Je )a r } 

^ R r ^3 ag g a op o op p 


(*-13) 
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Consider now : 


Jg 


or 


V R r e^ 4>,^ 

- J -a ka 

(R r ) 3 


-fr / (5 1 

(R > ii 


+ + r^ + 


— — rr f {a r * (R 1 + + r^ + q] $p)} {a r • i£} drn^ 

(R r ) 3 J a l l - 

B J 


1 ?« - - *-4 4 ... 4 ° £ - ^4 • 4 


(Rr) 3 ko oB 3 


(R r ) 3 



where 




R- 1 


+ 



Defining the following quantities: 



scalar 



one has 


v)8 


0 L _a i A je R r _ _L 

^ R r^3 ka a0 0 (R r ) 3 l " ka ” a & 


{«•?.. aII 


+ nj + mJ, qj - 3 A 3 ! a! <a r<T> A re(T) rr 1 ) 

k k£ ’A ka a0 0 0 crp p 


-3b 


J 


A J ! a! 


ka o0 0 


3 c kt ’i 1 


where 


1 . a r(T) Ja(I) B J 

ky o 06 k6y 


J = a r(T) A Je(T) j A je & r 
k£ o ca k£a0 0p p 


Finally, subtracting off the orbital term as described in Appendix C: 


v)G 


= - A 3 ' „ 3 + D 3 + M 3 q 3 - 3 4, q 3 

,„r s 3 ka a0 0 k Tt£ M £ k£ H £ 


on 


(B-14) 


- 3 a^ (a r(T) A*® (T) n j ) - 3 b:| A j * a£} 

ka a0 0 a op p ka a0 0 


Consider now Q <$ B : 


Q®i 8 


e jg 

-a 


/. 


{(r^ + u^) x d F^} 


= G' 


j( T > J . 

a0 -0 


/ 


{(r J + u^) x d F J } 
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But 


f {(. + u^) 


d p } = - S j6 X (R r + n j ) + (S j6 X a r )(a r • " n j ) 

(R r ) J (R ) 


+ 77^3 / [ ( ' J + “t ♦*> * i] [ ( ' J + "k *k> ' 2 r ] dmJ • 

(R > gj 1 1 


Using the relation (a x b) (c • d) - - b x (a c • d) , the last term above 

becomes 


- ^a r x 
<R C > 3 ' 


/ [? J + i J *i 4 + q i *i ;J + “1 +1 ♦£ , ’k] d " J • - r 


3 y m^ r 

= - ~T 3 ^ x 
(R ) 


a 3 + *l 4 + H m 4 + ^ e L% 3 )- 2 ' 


- - ^ UL ,- (- a r x I J£ • a r + a r x bj qj + a r x sj 11 qj + a r x S J „ 13 q J „> 


<R r > 3 


k q' + a x s k q k ■ r - - t a t - 


where: is the Body j inertia dyadic about its hinge 


gill , SJCT) . a t 


bJ „ A * 8 a 1 eJ 

Bo O ~a 


qj!3 _ j . qi qJ . ^ r 

S & ■ q k S kJl 5 S kfc C kJt - 


since ®^ 8 = 

Q®| 8 

- ' Y , . m A e J 

(R r ) 3 

• <S J6 x R r ) - - 

lusi e j . 

(S j6 x a r ) (a r 

(R r ) 3 -° 


3rai e J . 

a r x {- P f • 

(R r ) 3 ■“ 



(R r ) 3 “° 


r • Z 3 


(s ie x n^) 


n J ) 


J k ' “k “k ' M k 
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or 


®l* - - 


V -T. s j6 A je R r 

( R r ) 3 s «* A e* r * 


_ Y pA s j6 A jr n j 
(R r } 3 a3 A g6 n 6 


3_ymj cJ6 A je „r ,J(T) A re ^ 


(R ) 


r N 3 "ag "36 a 6 (n cr 


ap 


- A J * S aJ * (T) (- if A Ja a r + (bj™ + sf 1(T) + S J .J 3<T> )qJ}. 

, r.3 a3 3Y y6 yd a p p 6k 6k 6k k 

(R / 


Thus, subtracting off the orbital term as described in Appendix C: 


®f - 

- j 
3y m 

(R r ) 3 

A je 

«3 

~r 

a &Y 

A J ' (T) if A Je a r 
Y6 6 o op p 


- 

Y m 3 

(R r ) 3 

{ 5 J 5 

a3 

A jr 

36 

"1 - 3 ^ *1 <f (T) C ■? 

(B-15) 

+ 

3 

a3 

~r 

a 0Y 

A je(T) j(T) jll(T) s jl3(TK j } 

y6 ^ b 6k + b 6k b 6k ' V 



where 


s j6 ' = 

d^ + H 33 


a 

a a 


s jn _ 

ka 

r(T) je(T) 

a 6 A 6S 

s L 

s 3 12 = 

r (T) . Je (T) 

c J 

k£3 

a 6 A 6a 

kZa$ 

s j 13 = 

q j S J12 


&a 

q k kJla 
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B. 3 Derivation of Aerodynamic and Solar Pressure Produced Generalized Forces 


This analysis will be restricted to consideration of only a flat 
plate configuration. Thus, utilizing Equation (3) of Section 6.2, the 
general expression for the force on an element of area dA^ is given by 

dF^ F = P f |cos q^|{H^ F cos q 0 e 0F + 6 F } dA^ 

Therefore, the virtual work done by dP F acting on an arbitrary element 
dA^ due to the virtual displacement 

6R 1 + 6u^ + x (P + u j ) 


is 


where 


6W 


jf 


-/ 


dP f • [6R 1 + dP + 60^ x (P + P)] 


j 

R -i i r 

5R » 6R e 
a -a 


dP = 6q„ «|)o P 

l r la -a 


6 0^ 


50 j e 3g 

a -a 


Also, 


so that 


dp f = P F 6R 1 + Q 0p 6q 0 + QCE^ ? sP 

a a k k a a 


P F = f 

“ J A 


dF jf • e r 


-a 


4 F ‘ /. 


dP f • e J P 

-a Y ka 


l6l 


(B-16) 



q@i F - / 

,3 


dF^ • {e^ 8 x (r“* + u^)} 


-a 


jF iF rK IF 

where R , and Q Qp J are the generalized pressure forces associated 

"a ' i k a 


with the generalized coordinates R 1 , and 0 3 respectively. 


Thus , 


R 


J F - of 


= P^lcos IT 3 I {H^ f cos n 3 e 


Jf 


e r + G Jf 
-a 


e r } A J 
-a 


so that 


R jF = A j p f | cos n j| A jr(T) ^ R jf CQS n j B j 

ot otp p 


+ ^ 5 £ } 


(B-17) 


where 


Jf 


e o 

0 “0 


< 4 


In particular then, 


R^ A = pV 2 A J |cos n 3 |A^ (T) {2(l 


c^) cos Bg 


+ ° j t- it b 6*> * 


and with 


J = J £ . * a 


cos n = e 


then 


cos n J - B im (- A J ! bh 
a v a0 0 
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Therefore 


R' 


JA 


Also, 


- pV 2 A J A^ (I) {2 (1 - o 3 ) Eg (> 


+ o j A& bj) |B 3(I) A 3 6 b r | . 
36 6 1 a op p 1 


BjW A 3 * b r 
6 6y y 


R- 


JS - P s A 3 | cos n 3 |A 3 * (T) (2o 3 cos n 3 Eg 


+ (1 - V 3 ) Agg Sp 


where 6 D 

j 3 
cos n = 


i 6 6 0 5 

A„ S , so that with { S } = {6 } (the unit sun vector), 

3Y Y 1(T) ie e e 6 

B 31 J A 3 * S„ . Finally then, 

a a3 8 


jS = P S A j A Jr(T) '/oJ J fJ< T > aJ® c e 


a6 


(2v J B, 


( B e ( 


K s 

3Y Y 


) 


+ < 3 - ^ A « S 6> l B J <T) A op S pl • 


J* 


Now looking at 

<4 F - / *L <4 


P^|cos rr^|{H 3 ^ cos n 3 e 3 ^ + G 3 ^ 6^}dA 3 


J 


Since we are working with a flat plate, the following approximation is used: 

J 


so that 


/, 

A 3 

/ 


^ dA 3 - ^ 

ka j 

m 


/ +ko ^ 


B" 


dA 3 = A 3 0? 
ka 1 


ka 


(B-18) 


(B-19) 
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Then, 


4 F = ^ $ ka pf | cos cos e Jf + G^ f 5 f } 


and 


OA = - 


pV 2 A' 


3 4b {2<1 - ° 3 > B 6 


( bJ 5 (t> 4; »;) 


+ cf j aJ® b^}|Bj (T) A je b r | 
po o 1 a a p p 1 


and 


of - p s aJ 


«j s (2v 3 Bj (l 


J(T) ja a 

e py y 


+ - “ J > < s b , i b3(t> 4: 4i • 


Op p 


Now calculating q®? f , 


* f 


e 

-a 


is • rzi j . nJ 


(r J + U J ) X dF 




= e jg 
-a 


/ (r 2 + u j ) dA J x P f | cos rv^J { f cos 


+ G jf 6 f } 


dA^ (centroid of plate) , 


so that letting r -* = ~ / 

A^ J 

Q®f - A j ej g ‘ (r° j +qJ ^)x P f |cos n j |{H jf cos r, j e jf + G jf 6 f } 


= A^ P^ (r ^ H 3 ^ „ no r>3 4- p j f 1 I j I 

aa lr ay Y H C0S n + & y G } l cos H J | • 


(B-20) 


(B-21) 
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But, transforming from hinge axes to Body j axes. 


®r - Ms 


so that 


® jF = P f (r° j + H^ 3 ) {B j H^ f cos n 

v ay ay y 


F } | cos 


Thus , 


■ -» v2 aJ Kl + »i 6 3 ) ^ - “ 3 > b b ( b s 3<T) a s? ^ ) 


+ a 3 A^bJ}|BJ <T) A^b r | 

Bo 6 1 o a p p' 


and 


®f ’ pS a3 <^s + "«B )t2vJ B e( 


4 m a^ s b r 
o oy y 


+ (1 - v j ) Aj® Sj> | B J (T) A je b r | 

Bo 6 1 o op p‘ 


(B-22) 


(B-23) 
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Appendix C . Modification of UFSSP to Simulate The Dynamics 
of Rapidly Spinning Flexible Systems 

The material for this appendix has been taken from Reference 2 which 
contains a more complete description of the modification including test cases. 
The equations contained herein are presented in their original form with the 
important exception that only first order terms have been coded so that only 
the first order correction terms are presented in Equation (4.4) . 

1.0 INTRODUCTION 


The modification described herein extends the capabilities of the program to 
provide for the accurate simulation of rapidly spinning flexible spacecraft, 
flexible spacecraft having variable angular rates, and the structural 
dynamics of helicopter rotors. In general, the modification provides the 
capability of simulating the dynamics of rapidly spinning systems of bodies 
or systems with variable spin in the configuration of a topological tree 
having terminal flexible bodies. To use this special option of UFSSP, the 
terminal flexible bodies must be representible as space curves, however, 
the program will predict the general banding and torsional motion of the 
body including centrifugal stiffening effects. 


When a beam — a spacecraft appendage or helicopter rotor for example, 
rotates about an axis perpendicular to its own axis, the resulting centri- 
fugal forces have the effect of stiffening the blade to a certain degree 
thereby increasing the natural frequency of bending vibration of the beam. 
Traditionally, this effect has been accounted for by deriving a centrifugal 
potential energy term which is added to the elastic energy of the system 
or by summing the forces on portions of the beam mass. The present approach 
differs in that the appropriate terms are derived from a modified displace- 
ment function for the flexible body. That is, the original displacement 
‘ — * 

function for the body u is replaced by 


u 


-* 

u 


1 

2 



<£> • <#> 


dr(n) 

dn 


dn 


( 1 . 1 ) 


where r is the position vector of a point on the flexible space curve. 
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When used in the derivation of the kinetic energy, Equation (1.1), 
yields the approprate terms to account for centrifugal stiffening as well as 
for the effects of other components of acceleration. 

In the next section. Section 2, the geometry defining the terminal 
flexible space curves and their deflections is considered. In the following 
section. Section 3, the modifications to UFSSP equations of motion are derived. 

These results are summarized in Section 4 and the method of their incorporation 
in UFSSP is outlined. 

2.U GEOMETRY ASSOCIATED WITH DEFORMED SPACE CURVES 

In the present version of UFSSP, flexibility is modeled as a 
summation of separable functions of displacement and time. That is, the 
displacement of an arbitrary point r^ in Body j at time t is given by 

qj(t) P L (? j ) (2.D 

i=l 

The equations of motion as they presently exist include all terms 
which arise from spin and other components of system motion on the first 
order displacement, ^ (r) . Therefore, any additional terms to be included 
in the equations of motion must arise from higher order displacements 
Flexural shortening is the second order displacement of Body j and therefore 
logically must account for the additional terms required. This term arises 
because the displacement functions, , are defined over the entire surface 

of Body j and the distance along the deformed Gurve to a particular point is 
greater than the distance along the same path on the undeformed body. 

Consequently to insure inextensibility, the deformed shape u (r , t) must be 
corrected by an additional term. 

The correction term which is derived here is applicable to flexible 
bodies in the configuration of a spacecurve, that is, one dimensional flexible bodies 
More complex geometries, such as a closed loop or a surface, present additional 
difficulties in that the net distance traversed around every closed path 



167 



on the deflected body must be zero. This is a solvable problem whose solution 
mainly depends on the resolution of indexing problems and the establishment 
of suitable numerical algorithms. However, it is more complex in the 
case of structures of higher dimension than for space curves. 

However, restricting the modification to space curves yields results useful 
to many practical problems. 

Figure 2.1 shows Body j, a flexible space curve, with its 
undeformed shape defined by the vector function r(s), where s is the 
arc length along the undeformed body measured from the origin. The uncorrected 
displaced form of the body is defined by the vector function d(s,t) measured 
from the origin (a fixed point) along the deflected structure. The vectors 
r and d are related as follows : 

d (s) = r (s) + u (s,t) (2.2) 

— * 

where u (s,t) is defined in Equation (2.1) and the superscript has been 
d-ropped since it is unimportant to the discussion of this Section to distinguish 
Body j from other bodies. 



Body j , A Flexible Space Curve 
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Figure 2.1. 



— 1 

A small distance, As , along the deflected Body j the position 

vector is 


_ ! 
d (s + As ) 


r (s + As) 


_* 

+ u (s + 


As, t) 


( 2 . 3 ) 


An element of arc on the deformed curve is then given by 


As = d (s + As ) - d (s) 


( 2 . 4 ) 


or 


_ ' _ _ _* 

As = r(s 4- As) - r(s) + u (s + As,t) - u (s,t) 


( 2 . 5 ) 


and, in the limit 


_* 


ds 


' = d I(gI ds + AuJ.pt l dg 


ds ““ ds 

f 

The length of the arc differential ds is 


Idr 

dr + 

2 d? 



. du . 

— 

du 

du j 

I ds 

ds 

2 ds 

ds 

ds 

ds / 


ds 


( 2 . 6 ) 


( 2 . 7 ) 


where it is to be noted that 


£i - 1 


( 2 . 8 ) 


Since — is the unit vector tangent to the space curve at s and 
ds 

-* 

dr . du _ n 
ds ds 


( 2 . 9 ) 


that is, the incremental displacement is normal to the corresponding incre- 
mental length of the structure. 
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Using Equations (2.8) and (2.9), Equation (2.7) becomes 


The quantity 


ds = 

du t du 
ds ds 

ds 


I T* ~ 

Vi - £ 


ds 


ds ds 

«<<1, therefore we use the approximation 


(2 


' = A + 1 1*2. \ 

y 2 ds ' ds / 


ds 


( 2 . 


The differential flexural shortening is given by the following vector 
dr (s) 


(ds - ds ) 


ds 


- (-t 


du du , \ dr(s) 

77 ‘ 77 ds J “TT - 


_* 


( 2 . 


The original displacement function for the body, u (s), may 
be replaced by the following function 


— — df — JQ 

u (s,t) = u (s,t) + u (s,t) 


(2 


-* -* 


where 


-m , ' 1 f du’ du dr(n) , 

u (s.t) - - 2 J o dT ^ 


and n is a dummy arc length variable. 


Substitution of Equation (2.1) into the second term of Equation (2.13) 


yields 


n 


n 


"•E I <!<*> *lk 


(2 


i=l k=l 
s 


where 


ik 


1 f ( d$i djk \ dr , 

‘I 1 \W in 


10 ) 


11 ) 


12 ) 


.13) 


.14) 
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3.0 DERIVATION OF THE MODIFICATIONS REQUIRED TO THE UFSSP EQUATIONS OF MOTION 


As stated in Equation (2.13) the displacement function for the 

“ A 

flexible body, u (s,t) can be corrected for flexural shortening by adding 
a term u ^ (s,t) to u (s,t). The corrected displacement function is 

U j (s,t) = u*^(s,t) + u“^(s,t) (3.1) 

where - £ £ q | ^ 

i=l k=l 

This modified expression for u"^ may now be used to modify the 
expression for the position of an arbitrary mass point in Body j, p3 t 
and subsequently, the kinetic energy expression for Body j. 

The instantaneous position vector to an arbitrary point 

in Body j may now be written 

pj = R r + R* + t + r^ + u^ (3.2) 


where 


u J 


= u J 


-mj 
+ u J 


The velocity vector of an arbitrary field point in the flexible 
Body j is found from Equation (3.2) to be 


dp^ 

dt 


= R r + R + x R 1 + 

A • * « . * 

+ or x r** + u^ + up x u^ 


+ w 1 X l 


(3.3) 
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Some terms in Equation (3.3) are constant over Body j, it will 

simplify the process of evaluation of the kinetic energy integral to group 

these in one term. Redefine as 

dt 

= x^ + J x r^ + u^ x u^ 

at 


where 


it^ 1 3 R r + R 1 + i r x R 1 + + w 1 x 


(3.4) 


Substituting the expression for u^ , Equation (3.1), into 
Equation (3.3) we obtain 

-1 *1 

dp J dp J -mj -j -mj 

3i" ’ 3F + u + “ * " 


(3.5) 


where 


dp 




dt 


x^ 1 + o)^ x r^ + u ^ x u*^ , 


the previous definition of 


dp j 


dt 

The square of the velocity of an arbitrary point in Body j 


is then 


dai . 

dp j 

- . 

■ &* J + 2 #* J - 

dt 

dt 

dt 

+ 2 

iE* 3 

dt 

• (ur^ x 

G raj ) + ^ 

+ 2 

-mj 
u * 

(<i)^ X 

+ (u^ x u m ^) • (o)^ x 
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The last three terms In Equation (3.6) are fourth order in q , q and 
will thus lead to third order terms in the aquations of motion. Consequently, these 
terms may be neglected from this point. Using Equation (3.6), the kinetic 
energy of Body j may be written as follows: 


^ A 3 


1 f dp 3 dp J , j 

2 J dt ‘ dt dm 


6 


j 


- x* ♦ / if* 1 . dJ ♦ / g*. (5i , u"h dJ 

6 j 8 j 

Substituting the second of Equations (3.5) into Equation (3.7), 
we obtain an expression for the kinetic energy which is the sum of eight 
integrals : 

T » T* + x^ * J u mj dm 3 + x^ • y x J~ u™ 3 dm 3 ^ 


6 


j 


8 


j 


I 


@ 


& .J 1 )' dn»l + 


•/ = 


. G 

• J. J ' 


-*j -mj , j 
” J • u dm J 


/ (a ) 3 X u 


© . 


-j -*j. -mj , 

v .. J ) * u J dm 


^ + f (u> 3 x r 3 ) * (u) J x u'“ J ) dm J 


rZi „ aJ 


(3.7) 


f u*^ • (w^ x u m3 ) dm^ + f (m 3 x u ^) * (oj^ x u™ 3 ) dm 3 

6. 6, 


i 


i 


where T is the previous expression for the kinetic energy. 
From Equation (3.1) we have 


(3.8) 


-mj 

u 


E 4 < n 

i,k 


ik 


(3.9) 
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Consequently, 


f u mj dm^ = m - ^ 

J Q 


6 j 

” here 5 u)ik 

This quantity may 


q3 i ? U)ik 

i ,k 


m 


J" dm^ is a new special mass property, 

also be used to define the integral 


l ^ ^ - 2,i £ qJ i q k ?J d,ik 

J -I Ir 

1, tv 

We may now evaluate the kinetic energy term by term. 


tl = 2 m 


j -jl . V* J • j F j 

1 x L-*d q i q k (l)ik 


i,k 


t2 =■ m 


t3 = 


1 ;Jl ' (“ J x S q i q k f Cl>ik) 
i.k 

f (w 3 x r^) * U mj dm 3 

Jo 


-mj 

where Z., 
lk 


■ “ 3 • S "1 q k 1 

i.k 


ik 


- h f 

Jo 


m $ 


r 3 x F^, dm 3 
ik 


j 


The terms t4 and t5 will result in second order contributions 
to the equations of motion, consequently these are neglected. 


t6 = 


/ 

h 


(oo~^ x r^) ' (ur 5 x u raj ) dm 3 


m 3 -j 

T “ 


i.k 


4 4 ! ”k • “ 3 


( 3 . 10 ) 

(3.11) 

( 3 . 12 ) 

(3.13) 

(3.14) 


(3.15) 
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where 


■ h 1 < <fJ ' 5 kP ? - 5 «> d » j 


and 6 is the Kronecker delta. 

The remaining terms t7 and t8 are second order. 

The previous results may be combined to write the kinetic energy 
for Body j as follows: 


T = T 1 + T 2 


(3.16) 


where 


T 1 = T + a) 


" J ' (” J E q i q k ^ik) 


(E q i «£«)• 

' i.lc 


i ,k 


j J 


T 2 = m^ X jl • 


E 4 4 

i,k 


+ “ J x ^ q i 4 f («ik 


The terms in T.. are of the same form as certain terms in T , 

j. 


that is, in T we have the terms 


m -j 

— — U) 

2 


1 - (E q k ®ik) ' “ 3 + • (“ J E 4 4 H k ) 

i,k x i,k / 


(3.17) 


Thus, the terms in T^ can be accounted for in the original equations 
of motion by redefining and as follows: 
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(3.13) 


where 

and 


=*i 

E., 

xk 




eJ 

= E . + 




ik 

ik 

ik 



Z j 

-*1 

= Z J + 

z m , j 



ik 

ik 

ik 

and 

-*1 

Z ik 

replace 

E"). and 

ik 

zl as they are originally defined 

E ik 

9 

and 

ik 

, are modifications that must be added to 

rrect 

for 

flexural 

shortening 

• 


The terms contributed by to the equations of motion 

are established by means of Lagrange's Equations. 

4.0 SUMMARY OF FORMULATION RESULTS' AND APPROACH TO INCORPORATION IN UFSSP 


4. 1 Summary of the Equations of Motion 

The modified equations of motion are as follows: 


The q* Equation 


d /3T ' 

dt l aA 1 

\ 9T 

d / 3I i\ 8T i 

37 V 3qJ / 3qJ 


+ 

(- j £ 

A Jr R* - 2 FjJ 2 A^ 1 w 1 

ay . -y ka ay y a a 


+ 2 nr 5 F«P 

ka 

(A j6 R r + 

ay y 

*r>) 

(4.1) 

where f ^ 2 = 

ka 

J F j 

^e (l)Jfka 
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The R^ Equation 


d /3T \ ST = d / 3 M ST 

dt \SR i / SR 1 dt 'SR 1 ' SR"*" 


m 


— 

j A jr(T) j 3 

■ ae F g + 


2 f 4 2(t) a - 


gk 


F r 1 ^ 
By Y 


(4.2) 


+ 4 a) F^ 5 + U)| 

By y By Bp p 


] 


where F 


Jl _ 


a 




FtP = 2 F ^ 

Fd 2 <*£ q F (l)£ka 

j6 _ .j j2 

F a " q k F ka 


The 0^ Equation 


d / 

r ST \ 

ST 

d i 


S^ 

dt ( 

. *i / ” 

_„i 

dt \ 

l .„i/ 

„„i 

k3 V 

se x 


3e x 


+ m J F^ 1 (A^ 6 R r + A^ r R 1 
Aot cry Y ay Y 


A jl i 1 + S J7 ) 

aY Yd 0 a 


(4.3) 


The above relations will now be used to redefine the elements 
of the B matrix and C matrix defined by Equations (5-35) to (5-49) . 

* 

Old expressions for the sub-matrices will be designated by a superscript . 
(Note: only first-order correction terms are included here.) 
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*tl 3 - 2 L F i£ $ 
0=1 
3 

A14* . , 5 V' F j2 A jr 

B kg 2 2-f ku aB 

a=l 


r Ai , J J2 J7 

C k ' 2 “ F ko S ct 
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A31 

ct£ 


B^ 1 + 

al 


£ *« 


e=i 


j jrCf) F j2(T) 
oB W 


A32 „A32 


B = B 
ay 


ay 


,A33 _A33 


oB 


oB 


A34 A 34 


B ■ - B 

ag 


ag 


,a3 


= C 


A3 


(4.4) 


Various quantities used in Equation (4.4) are auxiliary quantities 
which are calculated at each time step. These quantities are generally 
expressible in terms of the special mass properties. In the following, the 
calculation of the special mass properties is clarified. 

4.2 New Special Mass Properties Required 

The UFSS Program has been modified to include additional terms 
for a terminal body which is flexible and in the configuration of a space 
curve. Any space curve may be approximated by a sequence of interconnected 
points, and consequently, the most logical representation of a space curve 
is in terms of a series of interconnected nodal points. All of the modifi- 
cations to the special mass properties are derived in terms of such a model. 
However, it may be desired to model a structure having significant lateral 
dimensions as shown in Figure 4.1 which requires several rows of nodes. 
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1 • » • • - • • ' - • 

6 11 16 21 26 31 36 


Figure 4.1. A Flexible Body Represented by 
Several Rows of Nodes 

The quantities associated with 'the rows of nodes defining a structure 
such as the structure of Figure 4.1 must be used to define the properties of 
the system represented by a single row of nodes as shown in Figure 4.2 




Figure 4.2. Nodal Numbering Scheme Admissable 
to the Analysis 
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The quantities associated with the structure of Figure A.l will 
be designated as follows: 


4m* J 

= the 

mass 

lumped at the £th node of Body j 

*j 

x a£ 

= the 

ath 

coordinate (a=l,2,3) of node 

Jt in Body j 

^aUk 

= the 

ath 

component of mode k at node 

l 

* 

= the 

total 

* 

number of nodes (nu < 64) 



These quantities are all input to the program in its original form, but the 
user is now constrained to number the nodes as shown in Figure 4.1. These 
quantities must be used to obtain equivalent properties for a space curve 
as shown in Figure 4.2. One additional input is required, the number of rows 
of nodes, however this can be accounted for by reinterpretation of the signi- 
ficance of the input flag designating whether a body is rigid or flexible. 

Thus, let f equal the number of rows of nodes, then the mass lumped at the 
t*l\ ^ 

i node in the space curve is given by 


W - £ 


Am 


*3 


£ = fj(i -1) + 1 


where ra. 


1 , . . . , m . 


* < 

= m . /f . 
3 3 


(4.5) 


Thus, for example, for the structure shown in Figure 4.2 we obtain 
the mass at the third node from the nodal masses of the structure in Figure 4.1 
as follows: 
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A»* J 

£ 

5 (3 - 1) + 1 


The coordinates of the nodes and the components of the modes of 
the representation shown in Figure A. 2 will be assumed to be those associated 
with the center row of nodes in the structure of Figure 4.1. Consequently, 
fj must always be an odd number. The number of the node in the representation 
of Figure 4.1, node £ , which corresponds to node i in the structure of 
Figure 4.2 is found from 

i * | (2 f i - f + 1) ( 

For example, for f^ = 5, and i = 3 we obtain 

£ = i (2(5)3 - 5 + 1) * 13 

The modes and coordinates of at node i in the representation shown 
in Figure 4.2 are found as follows: 

^oik a ^o£k * 


where £ 


i (2£ j 1 


fj + 1) and 


Preliminary to defining the modifications to the special mass 
properties, the following quantities will be required, none of which are 
output. 
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Ar 


ai 


x^. - x^.. a =1,2,3 

ai a, l-l 

i = 2,3, . . .m. 


A 1 j 

Ar J ., = x J , 
al al 


j _ 


( £ (Ar «i >2 ) 

a=l 


1/2 


and a vector is defined at each node, h , for modes k and l , 


ak£h 


z t <’ 

i -1 1 8=1 


j 

Ski 


“ ♦oi. t,_i) (♦]» - ♦ 


) 


8k,k-l VT eii y Bi,i-l / (4.9) 


Component B 
mode k 
node i 


h = 1, 2 • • • , m. 


The following quantities will be calculated. The first, 
is a special mass property which arises in base motion terms, while the second 
and third are terms which add to the existing special mass property terms. 
These additional special mass properties can be expressed in terms of 
integrals which take the same form as other existing special mass property 
integrals. 


The special mass property 


(l)k£ 


is defined as follows: 


^ / 5 k» dJ 

m R 

j 

This integral is replaced by a summation over the configuration of Figure 4.2, 
specifically by a summation over the i nodes defined by (-4.6). 
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Except for the dual subscripting, the special mass property, ’ 
defined in (4.10) is in exactly the same form as the integral defining 
in Equation (A-3), 


The second special mass property to be defined, 
modification to be added to an existing special mass property, 
is defined as follows: 

If £ 



a 


f 

v z 




P x <ta j 


This expression, except for the dual subscripting is in exactly 
the same form as the expression for the existing special mass property 
defined in Equation (A-4). 


The remaining special mass property, 
to be added to an existing special mass property, 
defined as follows: 


is a modification 


E J 


kZ 


E ia 


is 


E mj = — 

kZ J 
m 


/ ( 5 < ;J ' f kt) d " ;l 


B 


j 


This expression, except for dual subscripting, is in exactly the 
same form as the expression for the existing special mass property n| de- 
fined in Equation (A-7). 


(4.11) 


(4.12) 
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Appendix D. Modification of UFSSP to Calculate Dynamic Loads 

This appendix describes the methodology involved in calculating dynamic 
loads within the UFSS program. Determination of these loads can be divided 
into two basic phases. In the first phase, forces and moments acting at 
each interconnection between adjacent bodies are calculated within the UFSS 
program via information available from the original dynamics subroutines. 

In the second phase, the mode— acceleration method is utilized to calculate 
internal loads at any desired node point within a given terminal flexible 
body. These internal loads are calculated by a separate, stand-alone 
program operating upon a special loads history tape generated by the UFSS 
program. 

D. 1 Interconnection Force and Torque Calculations 

The synthesizing algorithm within UFSSP is based on elimination of the 
interaction forces and torques between adjacent bodies of the system model. 

By so doing, the final system of matrix differential equations includes 
only the unconstrained degrees of freedom of the spacecraft model. This 
retention of only the minimal number of degrees of freedom contributes 
greatly to reducing the cost of the solution for the dynamic response. 
However, when the dynamic loads are desired, the interaction forces and 
torques must be obtained explicitly. Specifically, they are calculated 
from the linear and angular velocities and accelerations and the coordinate 
transformations associated with the individual bodies of the spacecraft model 
by following the same sequencing algorithm used to synthesize the dynamic 
equations of motion. 


D.1.1 Interconnection Loads for a Flexible Body 

If Body j is a flexible body (with limb Body i) , its interconnection 
forces and torques can readily be determined using the equations from Appendix 
A. Specifically, from Equation (A-9) : 


F jl - A 1 ' 
aB 


a 


A3 , _A31 " r A 32 "j A33 

-c s + b bi q n + B Hr 9 y 65 


i tjA34 

“6 + B B6 R s 


(D-l) 
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Similarly, from Equation (A-12) : 


pJ 1 = 


= -C 


AZ 


+ B' 


A21 •• j A t,A 22 ~j 


al 


q~ + B 0' 
ay Y 


B af “8 + E tr B e‘ 


(D-2) 


The coefficient matrices in the above equations are available immediately fol- 
lowing the flexible load operation for Bodyj. The quantities q^ and are 
elements of the dynamic derivative vector depicted as {X} = {a} in Figure 8.1 
on^page 90; thus, they are available immediately following a call to the deriv- 
ative subroutine of the master program. The quantities and R* are not 
available for any bodies other than for Body 1; consequently, their explicit 
calculation must be added to the program. From Equation (A-5) : 



ji 

Ba 


• i 

0) 

a 


+ G 


j+ 

BY 


•‘1 

e J 

Y 


+ s 


J5 

6 


In a similar manner. Equation (A-32) can be rewritten to produce 


R j = R 1 + L l2 - L 32 i? 1 + 


ay 


'ay 


+ A lr<T > 

ay y ay Y 


where 


t J 1 = Jr(T) ~ j _ ; r A jr(T) 
ay aB By aB By 

-J r lJ 1 + A^ (T) y - 4% A 3 . r(T> . 
ay aB By ag By <*B By “B By 


(D-3) 


(D-4) 


(D-5) 


Thus , w 3 and R^ are computed recursively and can be considered as an addition 
’a a 

to the rigid auxiliary calculations. 

The above calculations for T 3 * and F^* are performed only at those time 

a a 

points called out for printing and/or plotting since the computations are only 
necessary for output. 


D.1.2 Interconnection Loads for a Rigid Body 

If Body j is a rigid body (with limb Body i),the interconnection forces 

and torques at its limb connection can also be determined using the equations 

from Appendix A. Specifically, from Equation (A-30) and the relationship 

W js = F^, immediately preceeding Equation (A-33): 
a a 


F^ 1 = -C A3 + B A31 
a a ax, 


«T*C 




r A33 rJ 
B aB V 


(D-6) 
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Similarly, from Equations (A-30) and (A-37): 


T^ — £ F^ 1 - C A2 + B A21 0*J S + B A22 + B A23 rJ . 

a aS 8 a aZ & a$ 8 ag B 


(D-7) 


The coefficient matrices in the above equations are available immediately fol- 

.. • A 

"1 s 

lowing the rigid load operation for Body j. The quantities 0 J are elements of 
the dynamic derivative vector, while and are as given by Equations (D-3) 
and (D-4) with the single change of replacing in (D-3j. 

As in the previous case of a flexible body, the above calculations for 
and F^are performed only at those time points called out for printing 
and/or plotting. 


D . 2 Internal Load Calculations 


In the UFSS program, flexible bodies are modeled in the traditional struc- 
tural dynamics sense as a system of joints (or nodes) which are interconnected 
by weightless finite element members (e.g., beams, plates...). All masses are 
lumped at the joints. The orthogonal functions used to describe the spatial 
deformation of the bodies are normally taken to be the orthonormal modes pro- 
duced by a standard structural dynamics program such as NASTRAN, SAMIS, or SMAP. 
In general, most such programs are based on small deflection theory, using the 
direct stiffness matrix finite element approach assuming linear stiffnesses. 

Such an approach allows for the generation of a specific transformation matrix, 
herein called the load transformation matrix (LTM) , by a systematic application 
of unit forces along each degree of freedom with all other forces set equal to 
zero. In particular, such a matrix allows internal member forces to be calcu- 
lated through the simple matrix multiplication operaton 


k, (t) = B F (t), 

&Q Q 


(D-8) 


In the above equation, each component of ( t ) represents a specific desired 

internal member load, B^q is the LTM and Fq ( t) is the applied force vector. 
Specifically, B^q contains, as rows, coefficients for each degree of freedom 
relating the desired load to a unit force applied successively along each degree 
of freedom of the structure (assumed to be three times the number of joints in 
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UFSSP since only the modal translations are used to describe deformations), 
and the ordering of the columns must be identical to the degree of freedom 
ordering within the mode shapes. The applied force vector, F^(t), contains 
both the externally applied forces (environmental forces, control forces...) 
and the inertial (d'Alembert) forces associated with the inertial velocity and 
acceleration of each joint. 

Since the LTM is generated as a static problem with no relationship to 
the joint, masses , the following development is intended to show how the iner- 
tial forces enter the calculations and exactly what model is to be used for 
generating the LTM. 

D.2.1 Basic Concepts 

Consider an arbitrarily moving flexible body, B, modeled in the tradi- 
tional structural dynamics fashion as a system of N discrete joints, each 
having lumped mass m 11 (n = 1, . . . , N) , with massless finite elements (beams, 
plates...) connecting the joints. As shown in Figure D.l, let e^ (a = 1,2,3) 
be unit vectors defining a right-handed, rectangular Cartesian reference frame 
fixed in the undeformed body (i.e., deformations of B are measured relative to 

to the e frame). Let the origin of e be denoted by 0, and let the external 

-*tx 

forces acting at 0 be denoted by F°. Let a 11 be the inertial acceleration of 
joint n, and let F n be the total external force acting on joint n. The gov- 
erning equations of motion for joint n are as follows: 



Figure D. 1 - Schematic of a Flexible Body B 
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K 

- .n -n p n \ q \ 

B a = F + 6 g , (D-9) 

k=l 

_ n J c 

where g is the force exerted on joint n by the kth structural element (here, 

_ 

K is the total number of structural elements and, obviously, g =0 if element k 
does not connect to joint n) . 

Summing over the entire body B: 

N * ( K . ) 

E n. n i“- E {?” + E 8 ( +5° 

n=l :-i ( k=l ) 

(note, if 0 is located at joint ri of B, then conventionally F 11 = 0, as its 
contribution is expressed in F°) 


But, 




K 

E 


k=l 




and, since the elements are massless, 

N 

y -nk 
^ g 
n=l 

Therefore, 

N 

E n -n 
m a 

n=l 

or 


0 . 


F n + F° 

n=l 



at each instant of time. 


(D-10) 


Consider the same flexible body, B, rigidly cantilevered at 0, as 
shown in Figure D.2. Let f n be the constant external force applied at joint 
n. Then, the static force equilibrium equations for B are as follows: 

N 

E f” + £° - o. 

n=l 

where f° is the reaction force vector acting on B at the cantilever point 0. 
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B 

• • • 



Figure D.2 - System B Cantilevered at 0 


If, at any instant of time t, one now sets 


< v n n ^ v 
f = F (t) - m a (t) r 


(D-ll) 


then, necessarily, 



and the internal loads within the cantilevered body shown in Figure D.2 are 
identical to those within the free body shown in Figure D.l at time t. 


One now generates the LTM for the system shown in Figure D.2 by success- 
ively applying a unit force along each translational degree of freedom with all 

other forces set equal to zero. Having obtained the LTM (B ), one determines 

up 

the desired internal loads within the system of Figure D.l at any time t as 
being equal to those present in the system of Figure D.2 when the force distri- 
bution (D-ll) is applied. Specifically, the internal loads L^(t) are determined 
as 


V'> ■ % F Q (t)> 


where F^(t) is a 3N-vector, with components 



(D-12) 
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D.2.2 Internal Loads Program 

Because of the core storage required to save the LTM and the fact that 
internal load calculations can be performed subsequent to the dynamic simulation 
of a given system, it is advisable that a separate progam "LpADS" be generated 
to calculate the internal loads within a given flexible body. Figure D.3 presents 
an overview of the program interfaces. As seen from Equation (D-12), only two 
basic quantities are required for the loads computation. The first quantity, 


is obtained from a structural dynamics program; the second quantity, F^(t), 
is obtained from the UFSS program via a special history tape. 


As previously mentioned, the row index, M> of B^q runs from one to A, 
where A is the total number of internal loads to be calculated; the column 
index, Q, runs from one to 3N, where N is the total number of joints in the 
structural dynamics model of the given flexible body, as input to UFSSP. Addi- 
tional description of the LTM is contained in the user manual for the L0AD 
program referenced under "Associated Documentation" on page of this document. 



Figure D.3. Basic Interfaces for Internal Load Calculations 
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Thus, given Equation (D-12), it only remains to prescribe the 
computation of F^(t); considering Equation (D-ll) , the problem is further 
reduced to determination of f n . Let us first determine a n (t). 

Differentiating Equation (A-2), one finds that for Body j: 

V n = J R r + R 1 + 2^u r x R 1 ^ + (u* x R 1 ^ 

+ u r x^oj r x R^ + x x (u 1 x 

+ l ^ + 2^ x 

+ 0 )^ x^r^ n +u^ n ^ ,+ u jn + 2to^ x u jn + x x^ n + u^ 11 ^ 
or, specifically identifying that part of a^ n which is joint independent. 

n = A Bj + x ( P n + u jn ^ + u jn 

+ 2oj^ x u jn 4- x [J x + 

= R* + 2^i») r x R^ +^u r x + oj r x^to r x R"^ 

+ m 1 + + m 1 x x + 2^w i x I 1 ^ 

In the above expression for A J , note the absence of the term R . As in 
the case of the derivation of the basic dynamic equations in Appendix A 
(see, for example, the discussion leading to Equation (A-7)...), the d'Alembert 
force acting on Body j , if Body j were located at the origin of the reference 
frame, has been removed from the acceleration term by equating it to the 
corresponding gravitational force as reflected later in Equation <D-15) . 

Finally, in component form, 


•s* 


where 




a Jn - + «J” + 2 wf 

a a a op B 


“aB 


♦ <) ♦ + «f) 


(D-13) 
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where 


w jn = j Jn 
W £ q Jl 9 £8 

_ "J 

B H £ 9 £8 

The external forces acting on joint n (F^ n ) will presently include gravita- 
tional forces and control thruster forces. (Geomagnetic, solar pressure and 
aerodynamic pressure forces as well as control torquers can be added at any 
time if a specific application calls for their inclusion.) From Equation 
(8-16), the control thruster force components at joint n are given by 


F> nC - A 3 ! Rj nT 
a aB 8 


(D-14) 


r i r j 

From Equation (B-ll) and the fact that a has components -A^ in the e^ 

frame and components (0 0 -1 } in the e^ frame, the gravity force components 

at joint n are given by 


p jnG = _ ynv 
“ (R 


— — T [a* + (r jn + W jn )l 

r ^3 |_ oiB 8 a a J 

|" R J + A J r( * ( r jn + \^ n ) 1 A j * 

(R r ) 3 L 3 3a * a J a3 


(D“15) 


m jn -r 

(Note absence of the term g ^ - which has been equated to the term 
m^ n R r arising from as privi^usly noted.) 

Finally, the external force acting on joint n is given by 


*J n = p jnC + F jnG 
a a a 


(D-16) 


and the total applied force components making up F^(t) For use in 
Equation (D-12) are given by 


- m jn a jn 


F - m- a“ 
a a 


(D-17) 


with F^ n and a^ n as given by Equations (D-16) and (D-13) respectively, 
a a 
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